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Abstract 


This  research  in  computational  electromagnetics  addressed  the  problem  of  predicting  the 
near-field  mutual  coupling  effects  in  phased  array  antennas.  It  developed  and  demonstrated  a  new 
analysis  technique  that  uses  the  finite  element  method  (FEM)  in  combination  with  integral 
equations.  Due  to  FEM’s  inherent  ability  to  model  inhomogeneous  dielectrics,  the  new  capability 
encompasses  many  radiator  types  that  were  not  amenable  to  analysis  by  previously-existing 
methods.  The  analysis  considers  the  general  case  of  a  radiator  in  an  infinite  array  that  is  fed 
through  a  ground  plane  by  one  of  three  types  of  waveguides:  rectangular;  circular;  or  circular 
coaxial.  Accurate  feed  modeling  is  accomplished  by  enforcing  continuity,  between  the  FEM 
solution  and  an  arbitrary  number  of  waveguide  modes,  across  the  ground  plane  aperture.  A 
periodic  integral  equation  is  imposed  at  a  plane  above  the  antenna’s  physical  structure  to  enforce 
the  radiation  condition  and  to  confine  the  analysis  to  a  single  array  unit  cell.  The  electric  field 
is  expanded  in  terms  of  vector  finite  elements,  and  Galerkin’s  method  is  used  to  write  the 
problem  as  a  matrix  equation.  The  Floquet  condition  is  imposed  as  a  transformation  of  the 
matrix,  which  is  equivalent  to  wrapping  opposing  unit  cell  side  walls  onto  each  other  with  a  phase 
shift  appropriate  to  the  scan  angle  and  lattice  spacings.  The  solution  of  the  linear  system, 
accomplished  using  the  conjugate  gradient  method,  gives  the  electric  field,  from  which  the  active 
reflection  coefficient  and  active  element  gain  are  calculated. 

The  theory  and  formulation  were  used  to  develop  a  general-purpose  computer  code.  The 
use  of  commercial  CAD  (computer-aided  design)  software  for  geometry  and  mesh  generation 
makes  the  code  geometry  independent.  It  was  validated  by  comparing  its  results  to  published  data 
for  arrays  of  open-ended  waveguides,  monopoles  and  microstrip  patches.  Predictions  for 
dielectric-clad  monopoles  were  validated  by  a  hardware  experiment.  Finally,  the  code  was  used 
to  predict  the  scanning  properties  of  arrays  of  printed  dipoles  and  printed  flared  notches. 
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PHASED  ARRAY  ANTENNA  ANALYSIS 
USING  HYBRID  FINITE  ELEMENT  METHODS 

I.  Introduction  and  Background 


1.1.  Introduction 

A  critical  problem  in  phased  array  antenna  design  is  that  of  controlling  the  mutual  cou¬ 
pling  between  individual  antennas,  or  radiators,  that  comprise  the  array.  Mutual  coupling  may 
reduce  antenna  efficiency  by  creating  a  reflection  mechanism  that  depends  on  both  the  radiator 
geometry  and  the  scan  angle.  The  limiting  case,  but  one  that  frequently  results  from  poor 
radiator  design  is  scan  blindness,  meaning  that  there  are  one  or  more  angles  in  the  desired  field 
of  view  where  the  reflection  is  total.  Since  mutual  coupling  is  inherently  a  near-field  electromag¬ 
netic  effect,  it  is  rarely  possible  to  achieve  an  acceptable  radiator  design  without  the  ability  to 
accurately  model  the  near  fields.  As  we  attempt  to  design  phased  array  antennas  for  new  applica¬ 
tions,  we  often  find  that  existing  field  computation  techniques,  most  notably  moment  methods, 
do  not  adequately  account  for  the  radiator’s  topology,  feed  structure,  or  dielectric  materials.  This 
is  especially  true  for  several  broadband  radiators  that  incorporate  dielectrics  either  as  structural 
support  or  as  electrical  loading  that  reduces  their  size.  Therefore,  improved  computation  tech¬ 
niques  are  required  for  use  as  design  tools  for  broadband  phased  array  antennas. 

The  objective  of  this  dissertation  research  was  to  develop  and  demonstrate  a  new  analysis 
method  versatile  enough  to  predict  the  performance  of  a  variety  of  radiators.  The  method  is  a 
hybrid  of  the  finite  element  method  (FEM)  with  integral  equation  continuity  conditions.  The 
integral  equation  for  one  of  three  types  of  waveguides  (a  summation  over  dominant  and  higher- 
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order  waveguide  modes)  provides  an  accurate  method  for  including  feed  structure  effects.  The 
use  of  a  periodic  integral  equation  to  represent  the  fields  above  the  array  enforces  the  radiation 
condition  and  makes  the  analysis  tractable  by  confining  it  to  a  single  array  unit  cell.  A  new 
periodic  boundary  condition  for  finite  elements  was  derived  to  account  for  the  mutual  coupling 
across  unit  cell  side  walls.  The  work  culminated  in  a  general-purpose  computer  code  that  suc¬ 
cessfully  predicted  the  scan-dependent  properties  of  a  variety  of  arrays  such  as  open-ended 
waveguides,  microstrip  patches  and  printed  flared  notches.  When  possible,  the  results  were 
confirmed  by  comparisons  to  data  in  the  scientific  literature  that  was  obtained  from  either  mea¬ 
surements  or  by  other  methods  of  calculation.  These  validation  cases  are  only  a  sampling  of  the 
capability  of  this  new  analysis  tool,  which  is  able  to  predict  the  scanning  performance  of  most 
array  radiators  that  are  in  use  or  proposed  for  use.  Its  future  use  as  a  design  tool  is  bound  to 
improve  the  performance  of  phased  array  antennas. 

1.2.  Phased  Array  Antenna  Electromagnetic  Analysis 

The  most  important  reason  for  near-field  electromagnetic  analysis  of  phased  arrays  is 
impedance  matching,  which  translates  directly  into  radiation  efficiency  and  low  VSWR  (voltage 
standing  wave  ratio)  to  prevent  damage  to  transmitter  components  in  high-power  applications. 
Paradoxically,  radiating  elements  that  are  well  matched  in  isolation  do  not  necessarily  remain  so 
when  they  are  arranged  in  a  closely-spaced  lattice  to  form  an  array.  Taking  one  radiator  to  be 
a  reference  element,  some  of  its  transmit  power  couples  directly  into  other  elements,  and  some 
power  from  each  of  the  other  elements  will  likewise  couple  directly  into  the  reference  element. 
The  power  returning  to  the  transmitter  by  way  of  mutual  coupling  represents  a  reflection  mecha¬ 
nism.  In  order  to  achieve  a  good  active  impedance  match,  each  element  must,  by  itself,  be 
slightly  mismatched,  so  that  its  self-reflection  vectorially  cancels  the  sum  of  couplings  from  other 
elements.  There  are  numerous  methods  for  controlling  these  mismatches,  such  as  altering  the 
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radiator  geometry  or  including  matching  circuits  in  the  feed  network,  but  they  can  only  be 
effective  if  an  accurate  solution  for  the  mutual  coupling  is  available  [1:1662]. 

Most  existing  solutions  for  array  mutual  coupling  use  the  method  of  moments  (MoM)  in 
conjunction  with  an  infinite  array  approximation.  The  approximation  makes  the  problem  compu¬ 
tationally  tractable  since  periodicity  conditions  may  be  used  to  restrict  the  analysis  to  the  space 
around  a  single  radiator,  called  a  unit  cell.  It  is  a  reasonable  approximation  for  large  arrays,  and 
it  is  also  used  as  part  of  the  present  method. 

1.3.  The  Need  for  Improved  Analysis  Methods 

There  are  several  trends  in  antenna  development  that  are  causing  array  element  designs 
tc  outstrip  the  available  analysis  methods.  One  is  the  desire  to  use  electronic  scanning  antennas 
in  applications  such  as  airborne  satellite  communications  and  airborne  surveillance  radar,  in  which 
the  antennas’  intrusion,  protrusion  and  weight  must  be  minimized  [2],  Another  is  the  growing 
trend  toward  millimeter  wave  frequencies,  leading  to  dimensions  so  small  that  it  is  impractical 
to  fabricate  and  assemble  and  array  one  element  at  a  time.  Yet  another  is  the  potential  cost 
reduction  of  using  printed  circuit  and  integrated  circuit,  or  monolithic,  fabrication  [3]. 

The  radiating  elements  that  are  proposed  to  meet  these  new  requirements  often  include 
irregularly-shaped  conducting  surfaces  in  combination  with  inhomogeneous  dielectrics.  An 
example,  shown  in  Figure  la,  is  the  "flared-notch"  element  [4],  The  dielectric  card  is  clad  with 
metal  on  the  back  side  except  for  a  slot  that  opens  progressively  wider  near  the  top.  The  front 
side  is  bare  except  for  a  microstrip  feed  line.  It  is  shown  here  fed  from  a  coaxial  transmission 
line  that  penetrates  a  ground  plane,  but  the  feed  line  could  be  microstrip  as  well.  A  recent  MoM 
analysis  of  this  radiating  element  in  the  phased  array  environment  used  the  simplification  shown 
in  Figure  lb:  the  dielectric  is  ignored;  and  the  feed  line  is  replaced  by  an  delta-gap  source  across 
the  slot  [5].  It  is  clear  that  this  model  cannot  predict  the  effects  of  high-dielectric-constant  sub- 
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(a)  (b) 


Figure  1 .  Flared  Notch  Radiator:  (a)  Printed  Circuit  Fed  from  Coaxial  Waveguide; 

(b)  Geometry  Model  for  Method  of  Moments 

strates  (monolithic  antennas  commonly  use  Gallium  Arsenide,  whose  relative  permittivity  is  12.8- 
12.9)  or  of  radiation  from  the  feed  structure. 

A  second  example,  shown  in  Figure  2a  is  a  dipole  radiator  metallized  on  one  side  of  a 
dielectric  card,  with  a  balun  feed  metallized  on  the  other  [6],  The  dielectric  may  be  trimmed  or 
notched  at  each  end  of  the  dipole  to  reduce  the  mutual  coupling  between  elements  located  at 
intervals  along  the  card.  Figure  2b  is  the  structure  actually  modeled  using  an  innovative  MoM 
approach.  In  this  case,  the  presence  of  the  dielectric  was  taken  into  account  by  using  a  Green’s 
function  for  a  parallel-plate  region  periodically  loaded  with  dielectric  slabs  [7], [8].  Some  results 
of  that  work  confirm  that  the  dielectric  has  a  pronounced  effect  on  the  wide  angle  scanning 
properties.  On  the  other  hand,  the  sweep-back  of  the  dipole  arms,  which  is  known  to  be  impor¬ 
tant  for  achieving  a  good  impedance  match  at  wide  scan  angles  [9).  is  neglected.  Also,  the  feed 
is  modeled  as  a  simple  delta-gap  source,  neglecting  the  effects  of  the  balun. 
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(a)  (b) 


Figure  2.  Printed  Dipole  Radiator:  (a)  Actual  Geometry  with  Microstrip 
Baiun  and  Coaxial  Feed;  (b)  Method  of  Moments  Model 


COAXIAL  FEED  GROUND  PLANE 


Figure  3.  Stacked  Patch  Radiator  with  Coaxial  Feed:  (a)  Two  Continuous 
Substrate  Layers;  (b)  Non-continuous  Top  Substrate 
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A  third  and  final  example  is  shown  in  Figure  3a.  It  is  essentially  a  rectangular  microstrip 
patch,  fed  through  the  ground  plane  from  a  coaxial  cable.  The  second  "proximity-coupled"  patch 
on  top  of  the  upper  dielectric  layer  is  intended  to  increase  the  overall  bandwidth  beyond  the  5% 
[10]  that  is  typical  of  a  single  patch.  The  array  properties  of  this  radiating  element  have  been 
predicted  using  MoM  [11],  with  the  coaxial  feed  represented  by  a  frill  current  source.  The  mo¬ 
ment  method  approach  for  this,  and  other  microstrip  problems,  relies  on  using  Green’s  functions 
for  layered,  infinite  dielectric  slabs.  Those  Green’s  functions  do  not  apply  to  situations  such  as 
Figure  3b,  in  which  one  or  both  substrates  are  not  continuous  layers. 

These  three  examples  illustrate  the  deficiencies  in  previous  mutual  coupling  computation 
techniques;  and  the  corresponding  new  capabilities  that  have  been  obtained  with  the  hybrid  finite 
element  method:  (1)  multiple  dielectrics  with  arbitrary  shape;  (2)  irregular  conducting  surfaces 
that  support  currents  flowing  in  arbitrary  directions;  and  (3)  detailed  feed  structures.  A  further 
objective  that  is  also  important  is  geometry  independence:  Whereas  each  of  the  three  examples 
discussed  above  used  a  specialization  of  MoM  to  the  particular  structure  and  required  develop¬ 
ment  of  a  separate  computer  code,  the  present  work  resulted  in  a  code  that  can  model  all  three, 
and  many  others  as  well. 

1.4.  Methods  in  Computational  Electromagnetics 

The  techniques  of  "classical  electromagnetics"  provide  formalisms  for  casting  physical 
problems  as  mathematical  boundary  value  problems.  Since  the  solutions  that  may  be  obtained 
by  the  purely  analytical  methods  are  restricted  to  canonical  geometries,  most  electromagnetic 
design  problems  of  current  interest  require  the  use  of  numerical  methods  to  obtain  a  solution. 
The  techniques  of  "computational  electromagnetics"  (CEM)  are  formalisms  for  mapping  boundary 
value  problems  from  continuous  to  discrete  forms  so  that  they  may  be  solved  by  computer. 
Therefore,  the  objective  of  CEM  is  to  produce  tools,  i.e.  computer  codes  with  which  device 
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designs  may  be  evaluated  without  resorting  to  hardware  experiments.  The  tools  must  have  four 
essential  characteristics:  effectiveness,  reliability,  efficiency,  and  versatility.  In  other  words,  they 
must  consistently  obtain  correct  results  at  reasonable  cost  for  a  variety  of  problem  geometries. 
The  previous  section  showed  that  the  shortcomings  of  current  methods  for  phased  array  near-field 
analysis  are  mainly  in  effectiveness  (due  to  simplifications  of  the  actual  problem  geometry)  and 
versatility  (due  to  the  restriction  of  each  to  very  specialized  geometries). 

The  requirement  for  effectiveness  limits  the  search  for  new  techniques  to  integral  equation 
(MoM)  and  partial  differential  equation  (PDE,  finite  element  and  finite  difference)  methods.  The 
former  is  by  far  the  most  well  advanced  for  solving  antenna  problems  because  the  integral 
equations  incorporate  Green’s  functions  that  satisfy  radiation  conditions,  forcing  all  valid  solutions 
to  decay  to  zero  with  increasing  distance  from  the  field  sources  (equivalent  currents).  The  PDE 
techniques,  on  the  other  hand,  more  easily  account  for  dielectric  inhomogeneities  for  which 
Green’s  functions  are  not  available.  The  finite  element  method  is  more  appropriate  for  devices 
with  irregular,  especially  curved  surfaces,  because  it  may  use  irregular  grids,  or  "meshes,”  while 
the  finite  difference  method  typically  uses  regular,  Cartesian  grids.  A  similar  need  to  model 
objects  that  include  inhomogeneous  dielectrics  has  led  researchers  in  electromagnetic  scattering 
to  consider  hybrids  of  the  finite  element  method  with  integral  equation  methods  [12]-[14],  By 
surrounding  the  object  with  an  imaginary  boundary  in  free  space  surrounding  the  scatterer,  the 
finite  element  method  may  be  used  to  solve  for  the  fields  inside  the  boundary  as  though  it  were 
an  enclosed  region,  and  an  integral  equation  is  imposed  to  ensure  field  continuity  across  the 
boundary.  Hence,  the  hybrid  finite  element  method  (HFEM)  appeared  to  be  a  likely  choice  for 
the  phased  array  radiator  problem  as  well,  provided  that  a  means  could  be  found  for  implement¬ 
ing  periodic  boundary  conditions.  The  success  of  that  implementation  and  a  demonstration  of  its 
benefits  are  some  of  the  important  results  of  this  dissertation  research. 
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The  next  chapter  will  pi  ^ent  an  overview  of  the  solution  approach,  which  involves  three 
novel  continuity  conditions  for  the  three  dimensional  finite  element  problem.  The  detailed  analy¬ 
ses  and  derivations  constituting  the  problem  "formulation”  (its  description  as  a  mathematical 
boundary  value  problem,  and  its  reduction  to  a  linear  system  of  equations )  are  given  in  Chapters 
III- VI.  They  are  intended  as  documentation,  or  as  a  trail  for  the  reader  who  would  attempt  a 
similar  solution  to  related  electromagnetic  problems.  They  are  not  essential  to  understanding  the 
results  of  validation  tests  and  hardware  experiments  presented  in  Chapters  VII-X. 

This  hybrid  finite  element  method  is  a  frequency-domain  approach.  Hence,  throughout 
this  document,  all  field  and  current  quantities  are  understood  to  be  time-harmonic,  with  the 
complex  exponential  time  dependence  suppressed. 
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II.  Solution  Overview 


2. 1 .  Problem  Description 

The  generic  problem  geometry  is  shown  (cross-section)  in  Figure  4.  The  "interior 
region,"  denoted  0,  is  a  section  of  an  array  unit  cell.  It  is  bounded  by  the  surface  denoted  T, 
whose  bottom  wall  is  the  ground  plane  plus  the  feed  waveguide  aperture.  The  top  wall,  called 
the  "radiation  boundary,"  is  an  imaginary  constant-z  surface  at  an  arbitrary  position  in  free  space 
above  the  radiator  structure.  The  side  walls  conform  to  the  unit  cell  boundaries.  For  example. 
Figure  5  shows  two  arrays  (stacked-patch  radiators),  one  with  a  rectangular  lattice  and  the  other 
with  a  triangular  or  "skewed"  lattice.  In  each  case  the  unit  cell  is  a  cylinder  extending 
indefinitely  in  the  +z  directions.  Its  side  boundaries  are  chosen  to  satisfy  the  periodicity 
conditions  (discussed  in  Chapter  VI  and  Appendix  C).  The  region  Q  is  formed  by  simply 
truncating  the  unit  cell  at  some  plane  above  the  array. 
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Figure  4.  General  Phased  Array  Radiator  Problem 


0  may  now  be  viewed  as  a  cavity  containing  any  number  of  material  regions,  each  with 
distinct  constitutive  parameters  e  and  n,  which  may  be  complex  (lossy).  This  research  will  only 
consider  linear  and  isotropic  materials.  There  may  also  be  voids  within  0  that  represent  the 
interiors  of  \  .rfectly  conducting  obstacles.  Infinitely  thin  wires  and  open  surfaces  such  as  patches 
and  strips  are  also  permitted. 

2.2.  The  Matrix  Equation 

The  solution  to  the  boundary  value  problem  represented  by  Figure  4  is  the  electric  field 
E(x,y,z)  everywhere  inside  Q  and  on  T.  HFEM  will  find  an  approximation  to  E  in  terms  of 
piecewise-continuous  expansion  functions  weighted  by  a  column  vector  of  coefficients,  denoted 
E.  Those  coefficients  are  the  sclut  ion  to  the  matrix  equation 

[SEE +  SEJ  }E  =  E™ 

Complete  explanations  of  the  terms  in  this  equation  are  given  in  succeeding  chapters,  but  briefly: 
The  matrix  SEE  is  sparse,  representing  local  interactions  between  field  sources  inside  0;  S^1 
represents  interactions  between  field  sources  on  the  nonconducting  parts  of  F  through  integral 
equations.  The  right  side  vector  Emc  is  due  to  a  field  incident  on  T  from  the  feed  waveguide. 
The  performance  parameters  that  are  of  greatest  interest  are  the  active  reflection  coefficient  and 
active  element  gain,  which  may  be  found  directly  from  those  parts  of  E  on  the  waveguide 
aperture  and  radiation  boundary,  respectively. 

2.3.  Finite  Elements 

The  region  0  will  be  subdivided  into  small  volume  elements.  Four-sided  tetrahedra  were 
chosen  because  they  conform  more  readily  to  irregular  and  curved  surfaces  than  other  popular 
choices  such  as  six-sided  "bricks."  The  volume  elements  are  often  referred  to  as  cells,  and  their 
four  vertices  as  nodes.  The  collection  of  tetrahedra  is  called  the  mesh.  Material  properties  will 


11 


be  assumed  constant  within  each  cell.  Figure  6  is  an  example  mesh  (one  quadrant  of  a  thick  disk) 
that  illustrates  an  important  flexibility  of  tetrahedron  meshes:  the  mesh  density  may  be  varied 
within  an  object.  Although  10-20  nodes  per  linear  wavelength  is  usually  an  adequate  sampling 
rate,  one  may  wish  to  sample  finer  in  regions  where  the  field  is  expected  to  have  singularities; 
or  in  order  to  capture  fine  details  of  an  object’s  geometry. 

Finite  elements  are  polynomial  functions  that  are  defined  over  individual  cells,  or  sub- 
domains.  Chapter  III  will  give  a  more  detailed  description  of  the  linear,  vector  finite  elements 
used  in  this  work.  These  functions  are  used  as  expansion  and  testing  functions  for  Galerkin’s 
method,  which  is  the  mechanism  used  to  reduce  the  boundary  value  problem  to  a  matrix  problem. 


2. 4.  The  Weak  Form  Functional 

The  principal  distinction  between  finite  element  and  moment  methods  (as  the  terms  are 
commonly  used  within  the  electromagnetic  research  community)  is  that  the  former  is  applied  to 
variational  statements,  while  the  latter  is  applied  to  integral  equations  ( 15:16].  The  variational 


Figure  6.  Subdivision  of  a  Volume  Region  into  Tetrahedra  (one  quadrant  of  a  disk): 
(a)  Interior  Edges  Visible;  (b)  Interior  Edges  Hidden 
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statement  used  here  is  the  weak  form  of  the  vector  wave  equation.  The  time-harmonic  form  of 
the  wave  equation  for  electric  fields  in  a  source-free,  inhomogeneous  region  is: 


VxJ_Vx£-^er£  =  0 

Mr 


(2) 


A  functional  is  constructed  by  taking  its  inner  product  with  a  trial  function,  W,  then  applying  a 
Green’s  identity  (see  Appendix  A): 


£(£)=[  _Lvx  W*  •  VxE-k\er  W*  •£ 


dv  +  jk0T]0  |  W*  -Jds  =  0  (3) 

r 


This  is  called  a  weak  form  because  the  Green’s  identity  has  shitted  one  derivative  from  the  field 
E  to  the  trial  function  W,  thus  weakening  the  differentiability  requirement  on  E.  This  functional 
has  three  difficulties  that  integral  equations  do  not: 

(a)  the  Helmholtz  equation  specifies  the  curl  only  and  not  the  divergence.  There¬ 
fore,  extra  effort  is  required  to  enforce  the  divergence  condition  V  •  (eE)=0. 

(b)  Boundary  conditions  are  not  included.  Thus,  although  (3)  is  not  restricted  to 
a  class  of  problems  with  uniform  boundary  conditions,  extra  effort  is  required  to 
ensure  the  satisfaction  of  all  boundary  conditions  that  are  present. 

(c)  The  radiation  condition  is  not  enforced. 

Chapter  III  will  discuss  how  (a)  and  (b)  are  resolved  by  choosing  vector  expansion  functions  that 
obey  the  divergence  condition  and  that  satisfy  boundary  conditions  at  both  pci  u.  ’  conductors  and 
dielectric  interfaces.  Chapters  IV  and  V  show  how  the  radiation  condm  ■  ;  forced  by  substi¬ 

tuting  an  integral  equation  for  J  into  the  boundary  integral  of  (3).  In  tli  >  c  of  the  waveguide 
aperture,  that  integral  equation  will  take  the  form  of  a  sum  over  waveguide  modes.  In  the  case 
of  the  radiation  boundary,  it  will  be  a  periodic  integral  equation,  which  may  be  written  as  a  sum 
over  spectral  domain  sample  points,  i.e.  Floquet  modes. 
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The  last  detail  to  be  addressed  is  the  problem  of  enforcing  periodicity  conditions  at  the 
unit  cell  walls.  The  implementation  is  straightforward  in  theoretical  terms.  It  has  been  accom¬ 
plished  by  others  for  two-dimensional  grating  problems  ( 16].  Chapter  VI  discusses  its  extension 
to  three  dimensions  and  the  means  for  implementing  it  algorithmically:  The  matrix  is  first 
constructed  as  though  the  unit  cell  walls  are  open  circuit  boundaries;  then  the  matrix  is  modified, 
creating  some  new  terms  and  removing  others.  The  effect  is  as  if  opposing  mesh  side  walls  were 
folded  around  onto  each  other  with  a  phase  shift  appropriate  for  the  array  scan  angle  and  the  unit 
cell  dimensions. 


2.5.  Development  Approach 

Figure  7  illustrates  a  generic  "cavity  array  problem,"  a  somewhat  simpler  problem  than 
the  "general  array"  problem  of  Figure  4.  It  is  still  a  phased  array  antenna,  but  the  radiators  are 
separated  from  each  other  by  conducting  walls.  Their  mutual  coupling  is  only  through  apertures 
in  a  conducting  ground  plane.  This  is  appropriate  to  a  restricted  class  of  radiators  such  as  open- 
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Figure  7.  A  Generic  Cavity  Array  Problem 


ended  waveguides,  horns  and  slots.  This  problem  embodies  all  the  same  aspects  as  the  general 
array  problem  except  for  the  periodicity  conditions  on  the  unit  cell  side  walls. 

Figure  8  is  a  further  simplification,  which  will  be  referred  to  as  the  "RF  device  problem." 
The  interior  region  is  again  a  cavity  with  perfectly  conducting  walls,  as  in  Figure  7.  But  here, 
both  inlet  and  outlet  apertures  lead  into  waveguides.  This  embodies  all  of  the  aspects  of  the 
cavity  array  problem  except  the  periodic  integral  equation. 

The  RF  device  problem  was  the  first  stage  of  this  research.  It  provided  validation  of  the 
approach  and  implementation  for  the  3D  finite  elements  and  the  waveguide  aperture  continuity 
conditions.  A  detailed  summary  of  that  work  is  given  in  a  separate  report  [17].  The  second 
research  stage  replaced  the  outlet  waveguide  modes  with  Floquet  modes  in  order  to  solve  the 
cavity  array  problem.  The  third  and  final  stage  in  the  algorithm  and  code  development  included 
the  periodicity  conditions  needed  for  the  general  array  problem.  Each  of  the  three  solutions  was 
validated  by  comparing  computations  to  results  published  by  other  authors,  obtained  by  methods 
other  than  FEM  (measurements,  mode  matching,  method  of  moments,  etc.). 
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Figure  8.  A  Generic  Passive  RF  Device  Problem 
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///.  Interior  Region  Problem  -  Finite  Element  Formulation 


The  first  part  of  the  problem  formulation  is  the  application  of  the  finite  element  method 
to  the  interior  region  Q.  This  will  ignore,  for  the  time  being,  the  field  continuity  conditions  on 
its  enclosing  boundary.  This  chapter  will  discuss  the  variational  form  of  the  problem  and  the 
restrictions  it  imposes  on  the  electric  field  approximation.  It  then  discusses  the  nature  of  the 
vector  finite  elements  and  shows  that  they  satisfy  the  restrictions.  Finally,  it  gives  the  derivation 
of  the  interior  matrix  terms  using  Galerkin’s  method,  and  their  reduction  to  algebraic  expressions 
using  coordinate  transformations  local  to  each  tetrahedron. 

3.1.  The  Variational  Statement 

The  boundary  value  problem  consists  of  the  operator  equation  (the  vector  wave  equation), 
boundary  conditions  and  applied  forces.  The  solution  is  a  function,  the  electric  field,  defined 
throughout  fl.  The  finite  element  method  attempts  to  solve  a  variational  equivalent  of  the 
problem. 

The  variational  statement  consists  of  a  functional  (usually  an  integral  containing  the 
unknown  function  in  the  integrand)  and  admissibility  restrictions  on  the  function  [18].  Admissible 
functions  are  those  that  are  in  the  domain  of  the  functional  and  satisfy  the  boundary  conditions. 
Appendix  A  discusses  the  two  forms  of  functionals  commonly  used  for  vector  electromagnetic 
problems  and  gives  the  rational  for  selecting  the  weak  form  (3). 

There  are  three  admissibility  restrictions.  First,  the  divergence  condition  V*(eE)=0 
is  necessary  to  ensure  a  unique  solution  since  the  operator  equation  only  specifies  the  curl  of  E. 
Second,  the  tangential  electric  field  must  vanish  at  the  surface  of  perfect  conductors,  i.e. 
AxE=0  .  Last,  at  interfaces  between  dielectrics,  tangential  E  is  continuous,  but  normal  E  is 
discontinuous,  i.e.  A*(e|E1)  =  ft*(e2E2)  .  It  will  be  shown  that  these  three  restrictions  are 
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satisfied  through  a  careful  choice  of  the  expansion  functions  used  to  approximate  E. 

3.2.  Scalar  vs.  Vector  Finite  Elements 

The  most  conventional  finite  elements  are  linear  functions  defined  relative  to  the  mesh 
nodes.  Within  a  single  tetrahedron  there  are  four  such  functions,  one  per  node.  Each  finite 
element  is  defined  within  a  single  tetrahedron  and  is  zero  everywhere  else.  Scalar  node-based 
expansion  functions  for  E  are  assembled  from  the  finite  elements.  For  example,  in  Figure  9, 
there  are  eight  tetrahedra  surrounding  the  central  node.  The  scalar  expansion  function  is  defined 
over  all  eight  cells,  is  equal  to  1  at  the  center  node,  and  goes  linearly  to  zero  at  all  surrounding 
nodes.  In  order  to  represent  a  vector  field,  one  choice  is  to  expand  it  in  terms  of  these  scalar 
functions  with  vector  coefficients: 

M 

s  =  1 


0.0 


.25 


< 1>  -  .50  K 

^  s 


.75 
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1.0  i 

'  7 

NODE  s 


•\  : 
-...  V 


-7. 


Figure  9.  A  Three-Dimensional  Linear  Expansion  Function:  (left)  Eight  Tetrahedra  Surroun¬ 
ding  a  Node;  (right)  Linear  Finite  Element  with  Surfaces  of  Constant  Function  Value 
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where  M  is  the  total  number  of  nodes  in  the  mesh.  However,  this  expansion  has  three  important 
disadvantages: 


(a)  The  boundary  conditions  at  perfect  conductors  are  difficult  to  enforce,  espe¬ 
cially  at  edges  and  tips  where  the  surface  normal  is  undefined. 

(b)  At  dielectric  interfaces  where  e  is  discontinuous,  the  electric  field  normal  to 
the  interface  is  discontinuous,  while  the  tangential  components  are  continuous. 

But  if  a  node  s  is  on  such  an  interface,  (4)  implies  that  aU  components  are  contin¬ 
uous.  Hence,  a  node-based  formulation  will  not  accurately  predict  the  field 
behavior  at  dielectric  boundaries  [19]. 

(c)  It  does  not  generally  satisfy  the  divergence  condition.  Failure  to  enforce  the 
condition  will  lead  to  spurious  non-physical  solutions  [20].  It  has  been  widely 
presumed  that  the  penalty  function  method  could  be  used,  but  Boyse  et.  al.  point 
out  that  penalty  methods  are  only  justified  for  positive  definite  functionals  [21], 
and  (3)  is  indefinite. 

Many  of  these  difficulties  can  be  circumvented  by  using  vector  finite  elements  in  an 
expansion  of  the  form 

N 

E  ~'EesVs(*>y>*)  (5) 

where  now  s  is  an  edge  index  and  N  is  the  number  of  edges  in  the  mesh.  The  particular  form 
of  sometimes  attributed  to  Nedelec  [22]  that  has  been  used  most  successfully  is  [23], [24] 

f.  -LitfMrfM)  (6) 

where  f;  and  fj  are  the  linear  scalar  finite  elements  defined  for  the  nodes  i  and  j  bounding  edge 
s.  Ljj  is  the  length  of  the  edge,  and  is  included  as  a  scaling  to  ensure  that  the  component  of 
tangential  to  the  edge  is  a  unit  vector.  Figure  10  illustrates  the  two  dimensional  version  of  this 
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Figure  10.  A  Two-Dimensional  Vector  Finite  Element:  (a)  Direction  (arrow  length)  indicates 
function  magnitude);  (b)  Constituent  Linear  Scalar  Finite  Elements 

vector  function  and  the  scalar  finite  elements  from  which  it  is  constructed.  This  choice  resolves 
all  three  of  the  difficulties  associated  with  the  node  based  formulation  since:  (a)  it  allows  a 
simple  yet  effective  means  for  imposing  conductor  boundary  conditions;  (b)  it  enforces  continuity 
of  tangential  field,  but  allows  the  normal  field  to  be  discontinuous  at  element  boundaries;  and  (c) 
it  is  free  of  divergence  since 


v  •  9,  =  vjr.  •  v/;.  v/t  -vfj-fjV-  v/. 

=o 


(7) 


(fj  and  fj  are  linear,  so  their  second  derivatives  are  zero). 

The  principal  disadvantages  of  the  vector  elements  are  algorithmic:  (a)  most  CAD 
software  generates  a  node  listing  that  must  now  be  converted  to  an  edge  listing;  (b)  the  direction 
of  each  edge  vector  must  be  accounted  for;  and  (c)  the  calculated  fields  must  be  converted  back 
to  vector  components  for  output  or  display.  One  further  disadvantage  that  has  been  cited  [19], 


19 


[23]  is  that  the  number  of  unknowns  is  larger  since  there  are  typically  4-5  times  as  many  edges 
as  nodes  in  a  tetrahedron  mesh.  Hence,  there  may  be  %  -  %  more  unknowns.  This  estimate  is 
excessive  for  three  reasons:  First,  it  assumes  three  unknowns  per  node,  but  four  are  actually 
necessary  to  achieve  an  effective  node-based  formulation,  using  vector  and  scalar  potentials  [25], 
Second,  there  are  no  unknowns  associated  with  edges  on  perfect  conductors,  since  the  tangential 
electric  field  must  be  zero,  whereas  all  three  components  of  the  electric  field  at  a  node  on  a 
conductor  could  be  nonzero.  Third,  the  node-based  formulation  may  require  finer  sampling  near 
conducting  edges  and  corners  to  compensate  for  the  uncertainty  in  the  direction  of  ft. 
Furthermore,  the  connectivity  between  edges  is  lower  than  for  nodes,  typically  by  a  factor  of  2, 
and  hence  the  number  of  matrix  entries  is  smaller  by  the  same  factor. 

Given  the  above  facts,  the  vector  finite  elements  are  clearly  the  better  choice.  The  next 
section  will  show  the  derivation  for  the  interior  matrix  terms  using  the  expansion  (5)  in 
conjunction  with  Galerkin’s  method. 

3.3.  Discretization  via  Galerkin’s  Method 

Galerkin’s  method  is  a  specialization  of  weighted  residuals,  in  which  the  trial  functions 
are  the  same  as  the  expansion  functions.  Its  use  is  permitted  when  the  expansion  functions  are 
in  the  admissible  space  of  both  the  direct  and  the  adjoint  problem,  as  discussed  in  Appendix  A. 

Substitution  of  the  series  expansion  for  electric  field  into  the  operator  equation  leaves  a 
residual  error  R  =  L(E)-L(E)  where  E  is  the  series  approximation  to  E  and  L  is  the  linear 
operator  [Vx  /ij'Vx  +  kger].  The  inner  product  of  R  with  a  trial  function  W  is 


f  R  -W*dv  =  F(E)  -  f  £  e, 
J  J  /=i 


V  X+i  ’  V  XW: 


k ler$t  •  W  *  jc/v  *jkQ  1,0  jy  •  w*ds 


This  includes  the  original  functional  since  <L(E),W)  =  F(E),  but  since  F(E)  =  0  from  (3). 


(8) 

The 
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weighted  residuals  procedure  forces  (R,W)=0  in  order  to  solve  for  the  coefficients  e,,  giving 

o  =  |  £  et  V-rV  v  xW*-^£rJ(  'W*  dv  +y*0T?0p-  iv*£/s  (9) 

0  ,=1  r 

Substituting  each  \£s,  one  at  a  time,  for  W  gives  N  equations: 

[  I  et  Pr  v  x$t  '  v  X^S_  *o  er^f  '^s  dv  +  jk0T)0  [j  '\frsds  =  0  (10) 

-J  /=!  J 

o  r 

The  order  of  summation  and  integration  in  (10)  n.ay  be  reversed  since  the  coefficients  e,  are  finite 
and  the  functions  and  VxJt  are  bounded.  Then  (1C  defines  a  system  of  N  equations  in  the 
N  unknowns  e,.  The  volume  integral  terms  are  the  entries  in  the  matrix  SEE  from  (1): 

sf/£  =  |  v  •  v  x#/-*o er$s  •$,  dv  (11) 

The  expansion  and  testing  functions  are  each  defined  only  on  the  collection  of  tetrahedra  adjacent 
to  the  corresponding  mesh  edges.  Hence  the  integration  is  over  Qsl,  the  collection  of  cells  shared 
by  edges  s  and  t.  These  matrix  equation  terms  will  be  computed  by  carrying  out  the  integrations 
in  (11)  analytically  using  a  transformation  to  homogeneous  coordinates. 

3.4.  Homogeneous  Coordinates 

The  homogeneous  or  simplex  coordinates  are  defined  locally  within  each  tetrahedron 
[26:266-274].  There  are  four  coordinates  t(,  t2,  t3  and  t4,  but  one  of  the  four  can  always  be 
eliminated  using  the  relationship  tj +t2+t3+t4=  1  .  The  coordinate  t;  of  a  point  anywhere  within 
the  cell  is  the  distance  to  node  i  from  the  opposing  face,  normalized  to  the  cell  height  along  that 
direction.  Hence  t;  is  equal  to  one  at  node  i;  and  zero  at  all  other  nodes  as  well  as  everywhere 
on  the  opposing  face.  The  transformation  is  given  in  terms  of  a  4x4  matrix  [T],  whose  elements 
are  the  16  cofactors  of  the  following  matrix,  U,  made  up  of  the  cell  vertex  coordinates: 
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u  = 


1  *1  >1 

1  *2  ?2 

1  *3  ?3 


*1 

*2 

^3 


(12) 


1  x4  y4  z4 

For  example,  T22  =  (y3Z4-y4Z3)  +  y,(z3-z4)+z1(y4-y3)  .  The  four  homogeneous  coordinates  are 
given  as  follows  in  terms  of  x,y,z: 


' 1 

1 

h 

_  [7-] 

X 

h 

6V 

y 

/4 

z 

(13) 


where  V  is  the  tetrahedron  volume.  These  coordinates  are  especially  convenient  since  the  scalar 
finite  elements  become  functions  of  one  coordinate  only: 


fi(x,y,z )  =  t{  =  ±[Tn*xTt2*yTn*zTn] 


Vfi  =  ±l*Tn+rra+tTu] 


(14) 


(15) 


and  the  limits  of  integration  are  simplified.  Most  terms  of  (11)  will  reduce  to  integrals  of  prod¬ 
ucts  of  two  scalar  functions,  which  have  the  simple  result: 


1  *_,1  *_,1  *2 

l\\fif)dxdydz  -  6 vjdt,  j  dh  |  ,,,jdi 3  -  |(1  *Siy) 


(16) 


cell 


where  5^  is  the  Kronecker  delta  and  i  and  j  may  take  on  any  values  between  1  and  4. 

3.5.  Volume  Integral  Computations 

The  volume  integral  computations  are  carried  out  by  visiting  each  cell  once  and  adding 
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its  contribution  to  Sf  f  for  every  pair  of  edges  that  are  part  of  the  cell,  excluding  those  edges  that 
are  on  perfect  conductors.  Let  i,j  and  m,n  be  the  node  indices  of  the  endpoints  of  edges  s  and 
t,  respectively,  1  ^i,j,m,n^4,  i?ej,  m^n.  s  and  t  are  global  indices,  but  i,j,m  and  n  are  local 
indices  defined  relative  to  a  cell .  Using  the  identity  V  x  (aVb) =aVx  Vb  + Vax  Vb  and  noting  that 
the  second  derivatives  of  the  linear  functions  fare  all  zero, 

VX?,  =  2  L'VfixVfj  (17) 

Vxf,  =  2  LtV/mxV/„  (18) 

Considering  the  first  term  of  (1 1)  separately,  note  that  the  gradient  terms  (17)  and  (18)  are  con¬ 
stants  and  may  be  taken  outside  the  integral.  Thus,  cell  k’s  contribution  to  the  first  volume 
integral  is 


yEE  1 


»(*) 


Mr 

MrtfV'*)4 


L'L'VfixVfj  •  V/mXV/n 


+  (^i'4^/2  “  Ti2Tj,){TmJn2  ~ 

+  (7i2^3-7i3^/2)^m27n3_7’m37n2)] 


(19) 


where  vk  is  the  volume  of  cell  k.  The  cell’s  contribution  to  the  second  volume  integral  is 


-fiWfr  vfm+fjfnvfi-  vfmyv 


Again,  the  gradient  terms  are  constants,  so  this  may  be  evaluated  using  (16).  The  result  is 


23 


(21) 


VEE2 

■W) 


*0  'rL,L, 

720V* 


E[(l+6iJTjirnl-(l+?>jJTilTnl 
-  ( 1  +  5in  )  Tj[Tml  +  ( 1  +  8jn )  TilTml  ] 


The  submatrix  SEE  in  (1)  is  the  sum  of  SEE1  and  SEE2.  Equations  (19)  and  (21)  are  a  closed- 
form  evaluation  of  the  volume  integral  (11),  written  as  an  algebraic  expression  in  terms  of  the 
geometry  of  a  cell  and  the  constitutive  parameters  contained  within  it. 

This  completes  the  discussion  of  the  interior  finite  element  solution.  The  next  two 
chapters  discuss  the  integral  equations  for  the  regions  exterior  to  I).  Those  integral  equations  will 
provide  expressions  for  J  in  terms  of  the  transverse  E  on  the  boundary. 
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IV.  Waveguide  Continuity  Conditions 


This  chapter  develops  the  integral  equations  for  the  RF  device  problem  illustrated  in 
Figure  8.  These  equations  apply  to  the  two  waveguides  entering  and  leaving  the  cavity  region, 
and  are  in  the  form  of  sums  over  waveguide  modes.  The  derivation  presented  here  is  generic, 
applying  to  any  type  of  waveguide  for  which  eigenmode  solutions  are  available.  Appendix  B 
gives  specialized  expressions  for  rectangular,  circular  and  circular  coaxial  waveguides.  Using 
the  mode  sum  form,  the  computation  of  matrix  entries  will  be  shown  to  reduce  to  inner  products 
of  the  mode  functions  with  vector  finite  elements. 

4. 1.  Combined-Source  Integral  Equation  and  Modal  Expansion 

The  derivation  of  the  integral  equation  uses  the  equivalence  model  shown  in  Figure  11. 
The  interior  problem  sees  zero  field  outside  the  cavity  region  and  equivalent  electric  and  magnetic 
currents  on  the  open  aperture.  Those  apertures  are  assumed  to  be  planar  and  located  in  the  end 
walls  of  the  two  waveguides.  They  do  not  necessarily  extend  across  the  entire  waveguide  cross- 
section  (irises  are  allowed).  The  exterior  problem  sees  zero  field  inside  the  cavity  and  oppositely- 
directed  equivalent  currents.  Notice  that  when  the  interior  and  exterior  problems  are  superposed, 
the  equivalent  currents  cancel  and  the  original  problem  is  recovered.  The  integral  equation 
applies  to  the  exterior  problem:  It  gives  the  fields  in  the  waveguides  in  terms  of  the  equivalent 
currents.  The  following  approach  is  similar  to  that  used  by  Harrington  &  Mautz  in  a  moment 
method  solution  for  open-ended  waveguide  radiation  [27J. 

Let  the  field  in  waveguide  A  (z<0)  be  comprised  of  a  unit-amplitude  incident  field  in  the 
dominant  mode  traveling  in  the  +z  direction,  plus  a  series  of  reflected  modes  traveling  in  the  -z 
direction.  The  transverse  (to  z)  electric  field  may  be  expressed  in  terms  of  a  complete  set  of 
orthonormal  mode  functions  g;  [28]: 
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(c) 

Figure  11.  Equivalence  Model  for  Waveguide/Cavity  Problem:  (a)  Original 
Problem;  (b)  Interior  Equivalent;  (c)  Exterior  Equivalent 


i=0 

The  subscript  t  denotes  transverse  and  the  7,’s  are  the  propagation  constants.  The  dominant  mode 
function  is  go-  The  sum  over  modes  i  includes  both  TE  and  TM,  as  well  as  both  sin(<£)  and 
cos(<£)  degeneracies  for  circular  and  circular  coaxial  waveguides.  The  complex  coefficients  C; 
are  unknowns  that  may  be  expressed  in  terms  of  the  solution  for  the  transverse  fields  in  the 
apertures,  e.g. 

ci  =  [  I  *  8ids  ~  50i  (23) 

\  z= o 

where  is  the  Kronecker  delta.  The  propagation  constant  for  mode  i  is  related  to  its  cutoff 
wavenumber,  kci,  by 
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(24) 


which  is  positive  imaginary  for  propagating  modes  and  positive  real  for  evanescent  modes.  The 
transverse  magnetic  field  is 


H?  =  Y0e'y°z(ixg0)  -  £  CiYiey‘z(ixgi) 

;=o 

where  \ ;  is  a  modal  admittance: 


Y, 


-jyt 

(TE) 

Jk0<r 

CTM) 

y^o 

I'1 

(TEM) 

The  boundary  condition  at  z=0  is  -JA  =  ftxHA  giving,  from  (25): 


-  -i-oso  *  E  c,y,g, 

*=0  «= 0 

Subsftuting  (23)  gives  the  final  form  of  the  integral  equation  for  waveguide  A: 


(25) 


(26) 


(27) 


E 

i=0 


[  %  |  ■ 

'  2=0 


2^o 


(28) 


Notice  that  the  equivalent  magnetic  current  is  involved  indirectly  by  EtA(z  =  0)  =  MAxz. 
Similarly,  the  integral  equation  for  waveguide  B  is 


E  tf*.'  (  £fl  ~JB  =  0  (29) 

The  primes  signify  the  fact  that  waveguide  B  may  be  a  different  type  or  size  than  waveguide  A. 
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so  it  may  have  different  mode  functions  and  modal  admittances.  The  right  hand  side  is  zero 
because  there  are  no  sources  in  wavegude  B. 

A  general  requirement  that  ensures  uniqueness  in  aperture  problems  is  that  continuity  of 
both  transverse  E  and  transverse  H  must  be  enforced  across  each  aperture.  This  usually  implies 
a  requirement  to  solve  for  both  transverse  field  components  independently  (or,  alternatively,  for 
both  M  and  J).  However,  the  nature  of  the  waveguide  mode  expansion  links  the  transverse  field 
components  through  modal  admittances,  that  is,  Ej  and  Ht  are  not  independent.  Therefore,  it  is 
only  necessary  to  solve  for  one  of  the  two,  and  the  obvious  choice  is  E,  for  consistency  with  the 
interior  solution. 

4.2.  Discretization 

The  integral  equations  (28)  and  (29)  give  expressions  for  J  that  are  substituted  into  the 
boundary  integral  of  (3).  Hence,  they  are  tested  using  the  same  trial  functions,  as  the  interior 
volume  integral  term.  For  compactness  of  notation,  let  ¥ f ;  and  'k®  denote  the  following  inner 
products: 

'8ids  (30) 

**  =  f  ^  'Z,ds  (31) 

Appendix  B  discusses  the  methods  by  which  these  integrals  are  computed  for  rectangular,  circular 
and  circular  coaxial  mode  functions.  The  testing  procedure  gives  the  equations 

eT  -2jk„0Y^%,seTA  (32) 
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(33) 


sf/=y'*oio£ ^ er, 

i=0 

sf,J  =  *Z,s,t€  rB  (34) 

i=  0 

In  practice,  the  mode  sums  may  be  truncated  at  32  or  less,  depending  on:  (a)  the  waveguide 
type;  (b)  the  nature  of  the  obstruction  at  and  near  the  aperture;  and  (c)  the  ratio  of  the  wave- 
number  to  the  cutoff  wavenumber  (more  modes  required  near  cutoff)  [17].  Note  that  there  is  an 
entry  for  every  pair  of  edges  s  and  t  that  share  the  same  aperture,  regardless  of  whether  or 
not  they  share  any  mesh  cells.  Hence,  this  matrix  is  not  sparse.  E‘nc  has  terms  for  all  edges  s 
that  are  in  aperture  A. 

4.3.  S  Parameters 

The  performance  of  passive  RF  devices  is  typically  expressed  in  terms  of  their 
"scattering,”  or  "S"  parameters.  They  may  be  found  from  those  coefficients  of  the  solution 
vector  E  that  correspond  to  mesh  edges  in  the  two  apertures.  The  modal  excitation  coefficients 
may  be  evaluated  from  these  as 

ci  =  E  es*?i~dOi  (35) 

c'i  =  E  e*  *!i  (36) 

,erB 

The  coefficient  C0  is  the  reflection  coefficient,  or  S1(.  The  transmission  coefficient  into  each 
mode  of  waveguide  B  is 
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When  there  is  only  one  propagating  mode  in  waveguide  B,  r0  =  S21  .  If  there  is  also  only  one 
propagating  mode  in  waveguide  A,  then  the  following  conservation  of  power  relationship  must 
hold:  |S„|2+|S21|2=1  . 

The  derivations  of  this  chapter,  combined  with  those  of  Chapter  III,  provide  a  framework 
for  a  computer  solution  for  two-port  RF  device  S  parameters.  The  computer  code  implementation 
and  validation  results  are  presented  later  in  Chapter  VII.  The  next  chapter  discusses  how  this 
methodology  is  extended  to  solve  for  the  properties  of  phased  arrays  of  cavity  radiators. 
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V.  Periodic  Radiation  Condition 


The  cavity  array  problem  (Figure  7)  was  formulated  as  a  straightforward  extension  of  the 
RF  device  problem  (Figure  8).  It  required  replacing  the  integral  equation  for  the  second  wave¬ 
guide  with  one  appropriate  to  a  radiating  aperture  in  an  infinite  array.  This  chapter  gives  the 
form  of  the  integral  equation  for  an  infinite  periodic  array  and  shows  how  it  is  reduced  to  matrix 
form. 

5.1.  Periodic  Integral  Equation 

The  equivalence  model  for  the  cavity  array  problem  is  essentially  the  same  as  Figure  1 1 . 
However,  the  aperture  on  the  right  side,  formerly  waveguide  B,  is  now  one  aperture  in  an  infinite 
uniform  lattice.  Therefore,  the  equivalent  currents  extend  indefinitely  over  the  outlet  aperture 
plane  in  both  x  and  y  directions. 

Each  radiator  in  the  array  is  assumed  to  be  excited  by  a  unit-amplitude  incident  field  in 
the  waveguide  from  z<0,  but  the  excitation  phase  may  be  different  for  each  element  in  order  to 
produce  a  beam  directed  towards  angles  60,  <t»0  in  spherical  coordinates.  The  phase  shift  as  a 
function  of  x  and  y  is 


*lx.y)  -  e-‘*-x  e'1*’’ 
\l/x  =  k  sin0o  cos0o 

ty  =  *  sin  00  sin  <£0 


(38) 

(39) 

(40) 


Figure  12  illustrates  the  notation  convention  for  an  array  lattice  with  an  arbitrary  skew  angle  y. 
The  aperture  shape  is  arbitrary.  The  fields  and  equivalent  currents  in  each  aperture  must  have 
the  same  magnitude  as  a  function  of  x  and  y.  A  mathematical  statement  of  the  phase  relationship 
in  (38)  is  Floquet’s  theorem  [29]: 
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y 


dycot7 


m  ■  -1  m  «  0  |  m  ■  1 

Figure  12.  Infinite  Array  of  Apertures,  Skewed  Lattice 


J{x+mdx+ndycoty  ,y+ndy)  =  J(x,y)e  i^mdx*ndycoiy'> e  J^yndy  (41) 

The  left  side  is  the  equivalent  current  in  any  unit  cell,  and  the  right  side  is  the  current  in  the  unit 
cell  centered  at  the  origin. 

Equivalent  currents  J  and  M  in  the  apertures  will  generate  vector  potentials  A  and  F  in 
the  half  space  z>d.  The  magnetic  field  due  to  these  is 

H(r )  =  V  XA  - ju  F  +  —  ‘f  (42) 

The  integral  equation  results  from  evaluating  Ht,  the  transverse  magnetic  field  in  the  plane  z=d 
and  using  the  boundary  condition  J  =  -2xHt.  Since  A(z=d)  is  entirely  z-directed  it  does  not 
contribute  to  Ht(z=d)  and  the  integral  equation  is 
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juz  x  F  +  —  VV  -F 


-7=0 


(43) 


F  is  an  integral  over  the  source  current  M,  and  the  limits  of  integration  are  infinite  in  x  and  y. 
Appendix  C  describes  how  it  may  be  transformed  into  the  following  infinite  summation: 


0  =  ~J(x,y)  +  £  £  Tmn-  EJkxmn,kymn)e 


~j^xmnX0  j^ymriy 


(44) 


m=-ocin  =  -ao 


k?  -k2  -K 

xmn  ymn 


■1/2 


2dxdykr) 


(k2-k2  )  k  k 

'  ymn'  xmn  ymn 

t  2 

^xmn^ymn  (*  ~^xmn^ 


(45) 


Euc  is  the  2D  Fourier  transform  of  the  transverse  unit  cell  aperture  field  and  kxmn  and  kymn  are 
sample  points  in  the  spatial  frequency  domain: 


2ir  m 


-k0slnd0  cos  <)>q 


(46) 


kymn  =  ^  -  lnm-CO±y  -  ko^e^o  (47) 

Uy  ux 

sometimes  referred  to  as  Floquet  harmonics.  The  summation  in  (44)  may  be  computed  numeri¬ 
cally  because  its  terms  decay  with  increasing  |m|  and  |n| ,  as  discussed  further  in  Section  5.3. 


5.2.  Discretization 

The  integral  equation  (44)  is  reduced  to  matrix  form  using  the  procedure  outlined  in  the 
previous  chapter:  First,  solve  for  J  and  substitute  it  into  the  boundary  term  of  (3);  second, 
substitute  the  series  expansion  for  E;  and  third,  substitute  each  \ps  in  turn  for  W.  Note  that  the 
integral  equation  involves  the  Fourier  transform  of  the  transverse  aperture  electric  field,  so  its 
expansion  will  be  in  terms  of  the  Fourier  transforms  of  the  ^,’s: 
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N 

£„  =  E 

/=  1 


(48) 


Zt(kx,ky)  =  JJ^|  ejkjX ejkyy dxdy 


(49) 


where  rR  denotes  the  radiating  aperture.  Testing  gives  the  matrix  terms: 


S.V-EE  (50) 

m  n  X 

1  R 

where  £tmn  denotes  it^xmn^mn)-  limits  of  integration  may  be  extended  to  ±  oo  since 
is  zero  outside  rR.  Then,  since  \j/s  is  a  real  function,  the  first  integral  may  be  recognized  as  the 
complex  conjugate  of  the  Fourier  transform  of  \j/s.  (The  Fourier  transform  of  a  real  function  is 
Hermitian,  i.e.  F(-k)=F*(k)  [30:193]).  Hence,  a  final  expression  for  the  elements  of  the  matrix 
is 

^St  =  %t(kXmn’kymn)  '  ^ mn  ’  %s  ^xmn,kymn)  ,  S,/  £  Ts  (51) 

m  n 


5.3.  Floquet  Mode  Limits 

The  infinite  sums  in  (51)  must  be  truncated  at  some  upper  and  lower  limits  ±m  and  ±n. 
Those  limits  are  easily  determined  from  the  form  of  £,  derived  in  section  C.4.  Figure  13  is  a 
contour  plot  of  | £j  in  dB  for  a  typical  finite  element.  The  axes  are  kxhx  and  kyhy,  where  hx  and 
hy  are  the  triangle  (mesh  cell)  heights  parallel  to  the  x  and  y  axes.  This  scaling  ensures  that  the 
size  of  the  contours  in  Figure  13  are  independent  of  mesh  cell  size.  The  Floquet  harmonics  are 
superimposed  as  dots  in  the  figure.  From  (46)  and  (47),  their  locations  (for  a  rectangular  lattice 
and  broadside  scan)  are  kxhx  =  2irmhx/dx  and  kyhy=27rmhx/dx  .  When  the  array  scans  away 


34 


15.00 


•SA 


Figure  13.  Typical  Scalar  Finite  Element  Fourier  Transform  (Contours  in  dB, 
sample  points  are  Floquet  modes  for  a  square  lattice  with  hx=hy=dx/5=dy/5) 


from  broadside,  the  points  will  move,  but  their  spacing  will  not  change. 

A  reasonable  upper  limit  on  the  number  of  sample  points  kxmn  and  kymn  that  must  be 
included  in  the  computation  of  (51)  are  those  inside  the  -20  dB  contour.  (The  product  of  £s  and 
£t  will  be  less  than  -40  dB  for  any  points  outside  that  contour.)  The  size  of  the  contour  is  consis¬ 
tently  kxhx  «  kyhy  *  ±2ir,  but  the  sample  spacing  is  inversely  proportional  to  the  unit  cell 
lengths  dx  and  dy.  Hence,  using  |kxhxj  ,  |kyhy|  <  2ir  in  the  leading  terms  of  (46)  and  (47) 
gives  ] m |  <  dx/hx  and  |nj  ^  dylhy  .  In  a  typical  problem,  dx  and  dy  are  each  approximately 
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.5X  and  the  mesh  ceils  are  approximately  .IX  in  height.  In  such  cases  limits  of  -5  <  m,n  < 
5  are  adequate  to  ensure  convergence.  More  modes  must  be  included  when  the  unit  cell  size  is 
larger,  or  when  the  mesh  cell  size  is  smaller. 

5.4.  Active  Reflection  Coefficient  and  Element  Radiation  Pattern 

Some  of  the  most  important  output  results  from  the  phased  array  analysis  are  the  active 
reflection  coefficient,  Ra,  and  the  active  element  radiation  (far  field)  pattern.  Both  are  functions 
of  scan  angle.  They  are  analogous  to  the  reflection  and  transmission  coefficients  in  the  two  port 
RF  device  problem. 

The  expression  for  Ra  is  identical  to  C0  in  (35).  Both  are  due  to  the  field  reflected  into 
the  inlet  waveguide,  computed  from  the  waveguide  aperture  field.  But  now  those  aperture  fields 
include  the  effects  of  a  periodic  radiation  condition  at  the  outlet  side,  and  so  it  is  the  reflection 
coefficient  for  one  feed  waveguide  in  an  infinite  array,  i.e.  the  "active  array  reflection  coeffi¬ 
cient.” 

The  active  element  pattern  is  analogous  to  a  transmission  coefficient.  It  is  a  measure  of 
the  excitation  strength  of  a  plane  wave  (a  Floquet  mode)  propagating  away  from  the  array. 
Amitay  et.  al.  show  that  the  6  and  0  polarization  components  of  the  element’s  far  field  pattern 
are  due  entirely  to  the  lowest  order  TE  and  TM  Floquet  modes,  respectively  [29:57).  Rewriting 
their  expressions  in  terms  of  Fourier  transforms  gives: 

Eg  =  (i cos#  +)> sin0)  -i*c(/cx0Q,ky00)  (52) 

y  dxdy 


E*  =  —=  .  (?cos4>  sin <t>)  ’E*c{kx00,ky00)  (53) 

V  dxdy 

A  check  for  conservation  of  power  may  be  made  by  computing  the  transmission  coeffi- 
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dents  for  the  two  mn=00  Floquet  modes,  whose  admittances  are  Y100  and  Y20o  f°r  TM  and  TE 
(see  section  C.5): 


Y0  is  the  feed  waveguide’s  dominant  mode  admittance.  When  the  feed  waveguide  supports  only 
one  mode  and  there  are  no  grating  lobes  in  visible  space,  the  conservation  of  power  relationship 
is  |Ra|  +  |Tfl|  +  |T*|  =  1  - 

The  derivations  given  in  this  chapter,  when  combined  with  the  preceding  chapters’  finite 
element  and  waveguide  derivations,  constitute  the  framework  for  a  computer  solution  for  the 
scanning  properties  of  cavity  arrays.  The  implementation  and  validation  results  are  discussed 
later  in  Chapter  VII. 
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VI.  Periodic  Boundary  Conditions 


The  final  stage  in  the  problem  formulation  bridges  the  gap  between  the  cavity  array 
problem  in  Figure  7  to  the  general  array  problem  in  Figure  4.  The  side  walls  of  the  cavity 
eg  ion  fi  are  no  longer  conductors,  so  adjacent  radiators  are  free  to  interact  across  those  walls. 
The  periodic  radiation  condition  does  not  account  for  that  interaction.  This  chapter  will  show  that 
the  general  array  problem  may  be  accomplished  by  constructing  a  matrix  for  a  single  unit  cell 
as  though  it  were  a  cavity  with  open-circuit  side  walls;  then  applying  a  mapping  directly  to  the 
matrix  to  enforce  the  periodicity  condition.  The  necessary  characteristics  for  a  unit  cell  will  be 
discussed  first,  then  the  algorithm  for  the  matrix  mapping  will  be  presented.  Next,  the  side  walls 
will  be  shown  to  have  no  net  contribution  to  the  boundary  functional.  The  specialization  of  the 
algorithm  to  edges  shared  by  the  radiation  boundary  condition  and  a  unit  cell  side  wall  is  dis¬ 
cussed  last. 

6. 1.  Unit  Cell  Representation 

A  typical  radiator  is  fed  by  a  waveguide  through  an  aperture  in  a  ground  plane.  Unlike 
the  cavity  radiators  considered  in  the  last  chapter,  it  has  some  structure  projecting  above  the 
ground  plane.  That  structure  may  be  enclosed  by  an  imaginary  box  whose  lower  surface  is  the 
ground  plane  at  z=0.  Its  top  surface  at  z=h  is  in  free  space  above  the  radiator  structure  (see 
Figure  4).  The  side  walls  of  0  are  the  unit  cell  boundaries.  As  was  indicated  in  Figure  5,  the 
unit  cell  may  be  trapezoidal  as  well  as  rectangular.  In  fact,  its  shape  is  fairly  arbitrary  within 
a  few  constraints,  with  some  possibilities  illustrated  in  Figure  14.  The  constraints  are:  (a)  the 
unit  cell  side  walls  do  not  cross  the  feed  waveguide;  and  (b)  opposing  boundaries  must  be  identi¬ 
cal  except  for  a  translation  of  (dx,0)  or  (dj,cot7,dy).  The  first  array  (Figure  14a)  has  circular 
waveguides  whose  diameter  is  larger  than  dy,  so  the  unit  cell  shape  along  the  boundaries  has  been 
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Figure  14.  Unit  Cell  Definitions:  (a)  Circular  Waveguide  Array;  (b)  Rectangular  Patch  Array 


altered.  The  second  is  a  rectangular  patch  array,  showing  that  the  unit  cell  walls  may  cut  through 
the  radiator’s  conducting  structure. 

The  unit  cell  definition  must  be  such  that  the  Floquet  condition  is  observed.  Let 
denote  a  field  in  the  m’th  column  and  n’th  row  of  the  lattice.  Then  the  unit  cell  fields  are  related 
by: 


_  /o 


,0 


.  „  ~JnnV 


ax  =  ^xdx 

ay  =  dy  C0t7  +  'Py  dy 


(56) 


where  and  \py  are  given  by  (39)  and  (40).  An  example  unit  cell  mesh  is  shown  in  Figure  15. 
The  two  perspective  views  show  opposite  sides  of  the  mesh  to  illustrate  the  important  requirement 
that  the  surface  mesh  on  opposing  faces  must  be  identical.  Every  edge  on  the  +x  or  +y  bound¬ 
ary  has  an  image  edge  on  the  -x  or  -y  boundary  respectively.  Consequently,  the  expansion  and 
testing  functions  associated  with  those  edges  are  identical  except  for  a  translation. 
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(a)  (b) 


Figure  15.  Example  Unit  Cell  Mesh  Showing  Opposite  Faces 

6.2.  Mapping  from  an  Infinite  System 

The  matrix  generated  by  discretizing  the  functional  represents  interactions  between 
electric  field  sources  associated  with  mesh  edges.  Consider  the  2D  triangle  mesh  shown  in 
Figure  16,  which  represents  an  infinite  mesh  of  an  infinite  array.  It  was  constructed  in  such  a 
way  that  opposing  unit  cell  boundaries  have  identical  mesh  edges.  Consequently,  when  the  mesh 
is  replicated  in  each  unit  cell,  the  finite  elements  are  perfectly  aligned  across  unit  cell  boundaries. 

If  the  finite  element  method  were  applied  to  the  infinite  problem,  it  would  generate  an 
infinite  size  matrix.  Appendix  D  shows  formally,  for  a  one-dimensional  problem,  how  periodic 
boundary  conditions  are  exploited  to  reduce  the  problem  to  a  finite  matrix  involving  a  single  unit 
cell.  The  extension  to  periodicity  in  two  dimensions  is  straightforward.  Essentially,  the  proce¬ 
dure  amounts  to  "folding"  opposing  unit  cell  boundaries  onto  each  other.  The  -l-x  and  +y 
boundary  edges  are  removed  from  the  problem  and  the  -x  and  -y  boundary  edges  will  interact 
with  those  just  inside  the  opposing  boundaries,  but  with  an  appropriate  phase  shift. 
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Figure  16.  Triangular  Finite  Element  Mesh  for  Infinite  Two-Dimensional  Periodic  Problem 


The  first  step  in  the  procedure  is  to  construct  the  matrix  for  an  isolated  unit  cell.  Next, 
the  terms  involving  the  +x  boundary  are  removed,  and  corresponding  -x  boundary  terms  are 
created.  Then  the  latter  step  is  repeated  for  the  +y  and  -y  boundaries.  For  example,  in  Figure 
16,  edges  15  and  19  share  a  mesh  cell  (a  triangle),  so  there  will  be  matrix  terms  S15  ,9  and 
S|9  |5.  Edge  2  is  edge  19’s  image,  so  the  new  entries  created  are  S15  2=S15  19exp{-jax}  and 
^2,i5=^i9,i5exP0O(x}  •  Also,  S19  19  will  be  added  to  S2  2-  This  is  necessary  because  of  the 
truncation  of  the  mesh  at  the  unit  cell  boundary.  In  the  infinite  mesh,  edge  2  would  interact  with 
itself  through  cells  to  the  left  and  right.  In  the  unit  cell  mesh,  there  is  no  cell  to  the  left,  so  S>o 
is  incomplete,  but  S19  19  is  identical  to  the  missing  contribution.  In  the  3D  problem  with  two- 
dimensional  periodicity,  there  will  also  be  matrix  entries  for  edge  pairs  that  are  both  on  the  +x 
or  +y  boundaries,  representing  their  interaction  within  a  cell  to  the  left  of  the  boundary.  These 
terms  must  also  be  added  to  their  -x  counterparts  without  a  phase  shift. 
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The  algorithm  is  summarized  in  Figure  17.  A  consequence  of  removing  mesh  edges  is 
that  some  of  the  matrix  rows  and  columns  become  entirely  zero,  and  they  must  be  deleted  before 
the  matrix  solution  can  begin.  This  is  consistent  with  the  fact  that  unknowns  associated  with 
image  edges  are  not  independent.  Due  to  periodicity,  once  one  is  known,  the  other  follows  from 
the  periodicity  condition.  Unless  one  set  of  dependent  unknowns  is  removed,  the  system  would 
be  over-determined. 

6.3.  Boundary  Functional 

The  algorithm  discussed  above  only  dealt  with  the  volume  integral  terms  from  the 

I.  FOR  EVERY  EDGE  sON  +x  BOUNDARY: 

A.  LOCATE  IMAGE  EDGE  s’  ON  -x  BOUNDARY 

B.  FOR  EVERY  EDGE  t  SUCH  THAT  Sst5*0: 

1.  IF  t  IS  ON  THE  +  x  BOUNDARY,  THEN: 

a.  LOCATE  IMAGE  EDGE  f 

b.  SET  Ss.,.=  Ss.,.  +  Ss( 

c.  SET  Sst  =  0 

2.  ELSE  IF  t  IS  NOT  ON  THE  +x  BOUNDARY,  THEN: 

a.  SET  Ss.,  =  Ssl  exp{jax} 

b.  SET  S|s,  =  Sts  exp{-jaj 

c.  SET  Ssl  =  Sls  =  0 

II.  REPEAT  I  FOR  +y  AND  -y  BOUNDARIES 

III.  COMPRESS  THE  MATRIX  (ELIMINATE  ZERO  ROWS  &  COLUMNS) 

IV.  COMPRESS  THE  INCIDENT  CURRENT  VECTOR 

Figure  17.  Periodic  Boundary  Condition  Algorithm 
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functional.  The  boundary  integral  must  now  account  for  the  contributions  from  the  unit  cell  side 
walls  in  addition  to  the  waveguide  aperture  and  radiation  boundary.  The  boundary  integral  may 
be  regarded  as  an  expression  for  power  flow  across  T.  In  an  infinite  array,  the  power  entering 
one  side  of  the  unit  cell  must  be  the  same  as  that  leaving  the  opposite  side,  so  it  is  expected  that 
the  side  walls  should  have  no  net  contribution.  The  following  will  show  this  to  be  true. 

Consider  the  boundary  integral  of  (3)  at  opposing  unit  cell  walls  parallel  to  the  x-z  plane. 
The  magnetic  field  must  obey  the  periodicity  condition,  and  since  the  outward  surface  normals 
are  opposite,  the  equivalent  currents  are 


J{x* 2/3,^)  =  -7( x,-^l)  e 

&  =  j  dy  coty 


-y«v 


(57) 


Any  admissible  trial  function  W  (see  Appendix  A)  must  obey  the  conjugate  relationship  because 
it  represents  waves  traveling  in  the  opposite  direction,  i.e. 

W*(jc  +  20,^)  =  ej0Cy  (58) 

The  boundary  functional  evaluated  at  the  two  unit  cell  walls  is 


y- 


h 

-]dz 


T-' 

f  W\x,-^)-J(x,-^)dx 

1-' 


(59) 


Fy+ 


*  \‘h  f 


W\x,^) -J{x,%)dx 


(60) 
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Using  a  change  of  variables,  let  t=x- 2/5: 


h  2  M 

Fy*  =  fdz  \  WlT+20,$)-J(T+20,$)dT  (61) 

0 

2 

Using  (59)  and  (60): 

dX  _o 

h  ~2  ^ 

Fy+  =  -\dz  |  W\t ,-d-l)*l{Ty-d-l)dT  (62) 

and  therefore,  Fy+  =  -Fy.  .  A  similar  procedure  will  show  that  Fx+  =  -Fx_  ,  hence  the  unit 
cell  side  walls  have  zero  net  contribution.  A  similar  result  has  been  reported  for  a  two-dimen¬ 
sional  problem  with  periodicity  in  one  dimension  [16]. 

6.4.  Radiation  Boundary 

A  form  of  periodicity  condition  was  already  imposed  on  the  radiation  boundary  through 
the  periodic  integral  equation.  However,  it  did  not  account  for  the  mesh  truncation  at  the  unit 
cell  boundary. 

In  order  to  complete  the  specification  of  periodicity  conditions  for  edges  on  rR,  either 
of  two  approaches  may  be  used:  The  first  method  is  to  apply  the  same  algorithm  (Figure  17) 
used  for  the  finite  element  boundary  terms,  but  excluding  I.B.2.  The  second,  demonstrated  by 
Gedney  for  periodicity  in  one  dimension  [16]  used  "overlap  basis  functions”  at  one  boundary 
edge.  Those  expansion  and  testing  functions  associated  with  points  on  one  boundary  extend  into 
the  next  unit  cell.  In  the  context  of  Figure  16,  that  amounts  to:  (a)  removing  the  expansion 
functions  for  edges  on  the  +x  and  +y  boundaries  (16-20);  and  (b)  extending  the  functions  for 
edges  1  &  2  and  3-5  into  the  next  unit  cells  to  the  left  and  below,  respectively.  Then  the  periodic 
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integral  equation  is  applied  to  the  modified  mesh.  Both  techniques  were  tested  as  part  of  the  code 
development  and  validation,  giving  equivalent  results. 

6.5.  Summary 

This  chapter  and  the  preceding  three  make  up  the  formulation  for  the  general  phased 
array  problem.  They  are  the  framework  for  mapping  the  boundary  value  problem  into  a  discrete 
form  that  may  be  solved  by  computer.  The  matrix  for  the  general  array  problem  is  constructed 
of  the  three  parts  that  were  developed  in  Chapters  III,  IV  and  V:  first,  a  sparse  matrix  due  to 
finite  element  interactions  within  the  unit  cell;  second,  a  dense  upper-left  submatrix  due  to 
waveguide  interactions;  and  third,  a  dense,  lower-right  submatrix  due  to  Floquet-mode  interac¬ 
tions.  The  side-wall  periodicity  conditions  are  implemented  as  a  transformation  of  that  matrix, 
as  was  described  in  this  chapter.  The  next  three  chapters  discuss  the  validation  tests  for  three 
computer  codes  that  implement  the  solutions  to  the  RF  device  problem;  the  cavity  array  problem; 
and  the  general  array  problem. 
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VII.  Validation  -  RF  Device  Problem 


The  goals  for  this  simplified  problem  (Figure  8)  were  to  demonstrate  essential  charac¬ 
teristics  of  the  finite  element  and  waveguide  mode  implementations.  The  3D  vector  finite  ele¬ 
ments  were  shown  to  correctly  predict  the  field  behavior  at  conductor  edges  and  dielectric  inter¬ 
faces.  The  higher-order  waveguide  mode  calculations  were  validated  through  comparisons  with 
measured  and  published  results  for  several  waveguide  discontinuity  problems.  These  results  also 
provided  estimates  for  sampling  requirements:  the  number  of  finite  elements;  and  the  number 
of  waveguide  modes.  More  complete  details  are  given  in  a  previous  report  [17]. 

7.1.  Computer  Code  Implementation 

7.1.1.  General  Procedure.  A  FORTRAN  computer  code  named  TWOPORT  implements 
the  solution  to  this  generic  problem.  An  outline  of  the  actions  it  takes  is  given  in  Figure  18.  The 
user  instructions,  read  during  the  first  step,  include:  the  type,  size,  location,  number  of  modes 
and  er  for  each  waveguide;  the  frequency  limits  and  frequency  stepsize;  and  the  name  of  the  file 
containing  the  problem  geometry.  The  geometry  file  contains  three  sections:  The  first  lists  the 
node  coordinates  and  several  flags  for  each,  identifying  boundary  nodes  (port  #  for  those  in 
waveguide  apertures  and  conductor  #  for  those  on  conducting  surfaces).  The  second  block  lists 
the  indices  of  the  four  nodes  comprising  each  tetrahedron  and  the  index  of  the  material  filling  it. 
The  last  block  lists  the  complex  er  and  ftr  of  each  material. 

The  third  step  converts  the  node-based  geometry  to  an  edge-based  geometry.  Each  edge 
in  the  mesh  defines  an  electric  field  vector  expansion  function  that  exists  over  all  cells  adjacent 
to  that  edge.  If  two  nodes  are  on  the  same  conductor,  then  there  cannot  be  a  field  along  the  line 
joining  them,  so  those  edges  are  not  included  in  the  edge  list.  This  is  the  means  of  enforcing  the 
boundary  condition  on  tangential  electric  field  at  conducting  surfaces.  On  the  other  hand,  if  two 
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I.  READ  INSTRUCTIONS  AND  OPTIONS 
n.  READ  PROBLEM  GEOMETRY 

III.  CREATE  EDGE-BASED  GEOMETRY 

IV.  FOR  EACH  FREQUENCY: 

A.  COMPUTE  TERMS  OF  SEE  ACCORDING  TO  (20),  (22) 

B.  FOR  WAVEGUIDE  A: 

1.  COMPUTE  INCIDENT  CURRENT  VECTOR  IN  (32) 

2.  COMPUTE  TERMS  OF  S^  FROM  (34),  (35) 

C.  FOR  WAVEGUIDE  B,  COMPUTE  TERMS  OF  S^  FROM  (34),  (35) 

D.  SOLVE  (SEE  +  Sn  )  E  =  Einc  FOR  E 

E.  COMPUTE  MODE  EXCITATION  COEFFICIENTS  FROM  (36),  (37) 

Figure  18.  Solution  Procedure  in  Program  TWOPORT 

nodes  are  on  different  conductors,  there  may  be  a  field  between  them,  and  such  edges  must  still 
be  included.  The  matrix  fill  operations  must  account  for  the  direction  of  the  vector  function,  so, 
by  convention,  it  is  always  directed  from  lowest  to  highest  node  index. 

Most  of  the  actions  under  step  IV  in  Figure  18  are  straightforward  implementations  of 
formulas  derived  in  Chapters  III  and  IV  and  Appendix  B.  The  matrix  elements  are  computed  in 
three  steps:  first,  the  interior  finite  element  interactions;  and  second  and  third,  the  exterior 
waveguide  interactions  for  waveguides  A  and  B,  respectively.  After  solving  the  system  for  the 
unknown  electric  field  coefficients,  excitation  coefficients  for  any  number  of  higher  order  modes 
may  be  computed. 

7.1.2.  Matrix  Solution.  The  matrix  structure  is  mostly  sparse,  except  for  two  dense 
submatrices  in  the  upper  left  and  lower  right  corners.  This  structure  results  from  ordering  the 
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edges  by  increasing  centroid  z  coordinate,  starting  with  those  in  the  waveguide  A  aperture  and 
ending  with  those  in  the  waveguide  B  aperture.  Unfortunately,  there  is  not  usually  any  special 
structure  (sparsity  pattern)  to  the  matrix.  It  tends  to  be  banded,  but  the  bandwidth  is  large  and 
there  is  no  advantage  to  using  band  storage  or  specialized  band  solvers.  Hence  the  two  approach¬ 
es  used  for  the  matrix  solution  were:  (a)  ordinary  LU  decomposition  (LUD)  using  standard 
library  routines  and  ordinary  row-column  storage  (IMSL  [31:34]  and  LAPACK[32:150])  for 
problems  with  2000  unknowns  or  less;  and  (b)  the  conjugate  gradient  method  (CGM)  based  on 
formulas  from  Sarkar  and  Arvas  [33],  using  sparse  storage,  for  larger  problems. 

The  CGM  solver  was  written  specifically  to  accommodate  the  form  of  sparse  storage  most 
attractive  for  this  particular  class  of  problems.  The  two  dense  submatrices  due  to  the  waveguides 
are  stored  in  ordinary  row,  column  format;  while  the  sparse  finite  element  matrix  is  stored  in  a 
column  array.  Its  entries  are  in  arbitrary  order,  stored  in  the  order  that  each  edge  pair  is  first 
encountered  as  the  fill  algorithm  performs  operations  one  cell  at  a  time.  That  assembly  technique 
is  considerably  more  efficient  than  performing  the  same  operations  one  edge  at  a  time  [26]. 

Figure  19  is  a  comparison  of  execution  times  versus  number  of  edges  for  LUD  and 
CGM.  The  CGM  solver  clearly  has  the  advantage  for  problems  with  1000  unknowns  or  more. 
Unfortunately,  its  solution  time  is  more  difficult  to  predict  a  priori,  since  the  number  of  iterations 
it  will  need  to  converge  is  unknown.  The  convergence  measure  is  the  residual  error  norm,  which 
is  the  L2  norm  of  the  solution  error: 

e,  =  |  SEi~E *"c  ||  (63) 

(i  is  the  iteration  number).  The  initial  guess  E0  is  the  zero  vector.  Figure  20  is  an  example 
convergence  history  (for  a  microstrip  transmission  line  fed  by  coaxial  waveguide  at  each  end). 
The  S  parameters  have  converged  to  within  1  %  of  their  correct  values  when  eje0  is  less  than 
.001 .  Problems  with  rectangular  and  circular  waveguide  ports  typically  only  require  1 N-2N  itera- 
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Figure  19.  Comparison  of  Solution  Times  for  LU  Decomposition  (IMSL)  and  Conjugate 
Gradient  Method  (CPU  seconds  on  VAX®  8650  minicomputer) 

tions  to  converge,  where  N  is  the  number  of  unknowns. 


7.2.  Waveguide  Discontinuities 

7.2.1.  Iris  in  Coaxial  Waveguide.  Two  test  cases  are  representative  of  the  several 
waveguide  discontinuity  problems  discussed  in  [17]:  a  conducting  iris  in  coaxial  waveguide;  and 
a  step  discontinuity  in  circular  waveguide.  For  the  first  of  these.  Figure  21  is  the  tetrahedron 
mesh  used  as  the  input  to  program  TWOPORT.  Relating  this  model  to  Figure  8,  the  inlet  and 
outlet  waveguides  are  both  coaxial  (inner  radius  a=  1.5mm,  outer  radius  b=3.5mm)  and  the 
cavity  region  fl  represented  by  the  mesh  is  simply  a  short  section  of  the  same  waveguide. 
Shading  has  been  added  to  identify  the  nodes  tagged  as  perfect  conductors.  One  quadrant  of  the 
inlet  end  is  blocked  by  a  thin  conducting  iris. 

Measurements  of  this  "device"  were  made  by  inserting  a  foil  iris  between  two  APC-7mm 
adapters  and  using  a  network  analyzer  to  obtain  S j }  over  a  2-18  GHz  frequency  range.  Figure 
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RESIDUAL  NORM 


Figure  20.  Convergence  of  Residual  Norm  and  Reflection  Coefficient 
using  Conjugate  Gradient  Matrix  Solver 


Figure  21.  Finite  Element  Mesh  for  Coaxial  Waveguide  Section  with  Conducting  Iris 
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REFLECTION  COEFFICIENT 


22  is  a  comparison  of  those  measurements  with  TWOPORT  calculations,  for  both  the  magnitude 
and  phase.  These  results  validate  two  important  features  of  the  solution  approach:  First,  the 
higher-order  coaxial  mode  functions  and  their  inner  products  with  finite  elements  are  correctly 
implemented.  Second,  the  finite  element  solution  is  correctly  predicting  the  field  behavior  at  the 
edge  of  a  perfect  conductor. 

7.2.2.  Circular  Waveguide  Mode  Converter.  A  step  discontinuity  in  a  circular  wave¬ 
guide  is  a  simple  type  of  mode  converter  commonly  used  in  feeds  for  reflector  antennas  [34], 
The  antenna’s  radiation  pattern  shape  is  controlled  by  carefully  adjusting  the  amplitude  ratio  of 
modes  in  an  oversize  (multimode)  waveguide  or  horn.  Figure  23  shows  finite  element  meshes 
for  two  test  cases  with  different  inlet  waveguide  radius  and  the  same  outlet  waveguide  radius. 
Multimode  calculations  for  these  geometries  were  presented  by  Masterman  &  Clarricoats  [35] . 
Their  results  are  in  terms  of  a  mode  conversion  ratio,  M: 


2  4  6  8  10  12  14  16  18 

FREQUENCY  (GHz) 

Figure  22.  Measured  and  Computed  Sn  Magnitude  and  Phase  for  Coaxial  Iris 
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M  = 


(64) 


I  (jgp2(P=a^  I 

I  c\s'px^P  =  a)  ! 

C|  and  C\  are  the  excitation  coefficients  for  the  TEt  1  and  TMj ,  modes,  respectively,  and  g'pi  and 
g'p2  are  the  radial  components  of  the  mode  functions.  Hence,  M  is  the  relative  strength  of  the 
two  modes  measured  at  the  wall  of  the  outlet  waveguide. 

Figure  24  compares  TWOPORT  calculations  with  the  multimode  results.  The  disconti¬ 
nuity  at  6.25  GHz  for  aj  =  1 . 15"  is  due  to  the  fact  that  the  inlet  waveguide  also  supports  the  TMj , 
mode  above  that  frequency.  The  agreement  of  the  TWOPORT  calculations  with  the  well-estab¬ 
lished  multimode  calculations  indicates  that  the  higher-order  circular  waveguide  mode  functions 
and  their  finite  element  inner  products  have  been  implemented  correctly.  It  also  demonstrates 
a  flexibility  of  the  approach  of  using  higher-order  modes  in  modeling  feed  structures  that  support 
more  than  one  propagating  mode. 

7.3.  Printed  Circuit  Devices 

One  of  the  most  important  capabilities  of  the  finite  element  method  is  its  ability  to  deal 
with  inhomogeneous  dielectrics.  Printed  circuit  devices  are  an  example  application  that  is  espe¬ 
cially  interesting  to  this  research  because  their  modeling  requirements  are  similar  to  the  printed 
antennas  discussed  in  the  introduction.  The  specific  test  cases  were:  (a)  a  straight  length  of 
microstripline  with  coaxial  connectors;  and  (b)  a  microstrip  meander  line. 

7.3.1.  Microstrip  Transmission  Line.  The  finite  element  mesh  shown  in  Figure  25  is  a 
thin  dielectric  slab  (100  pm),  situated  in  the  bottom  of  a  perfectly  conducting  box.  The  mesh  for 
the  air  space  above  the  slab  is  not  shown.  Shading  has  been  added  to  identify  the  nodes  associat¬ 
ed  with  microstrip  line,  the  coax  center  conductor,  and  the  coax  dielectric.  There  is  an  identical 
coaxial  port  at  the  far  end  of  the  device.  The  coax  dimensions  (a=43  /on,  b=  100  pm)  were 
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MODE  CONVERSION  RATIO  (dB) 


Figure  23.  Finite  Element  Meshes  for  Circular  Waveguide  Step  Discontinuities 
(Mode  Converters)  with  a2  =  1.4":  (a)  a,  =  1.05";  (b)  a,  =  1.15". 


FREQUENCY  (GHz) 


Figure  24.  Mode  Conversion  Ratio  for  Circular  Waveguide  Mode  Converters: 
TWOPORT  and  Multimode  [35J  Calculations 


CASE  OUTLINE 


Figure  25.  Microstrip  Transmission  Line  Section  in  Conducting  Enclosure 


chosen  strictly  for  convenience  since  the  objective  of  this  test  was  to  verity  a  microstrip 
transmission  line  property  (guide  wavelength). 

The  substrate  dielectric  constant  was  chosen  as  12.9  to  represent  Gallium  Arsenide 
(GaAs).  The  75  /on- wide  microstrip  line  has  the  same  characteristic  50  0  impedance  as  the  coax. 
The  enclosure  dimensions  were  chosen  to  ensure  that  there  are  at  least  4  line-widths  space 
between  the  microstrip  and  the  cavity  side  and  top  walls.  That  is  adequate  to  ensure  that  the 
enclosure  does  not  influence  the  guide  wavelength  or  characteristic  impedance  [36], 

The  transmission  coefficient  for  a  line  length  of  500  fun  was  computed  first.  Next,  the 
calculation  was  repeated  with  finite  element  geometry  scaled  by  1.5,  giving  a  line  length  of  750 
/xm.  The  difference  in  S21  phase  was  35°  at  the  40  GHz  test  frequency,  from  which  the  guide 
wavelength  is  calculated  as  2.57  mm.  Quasi-static  formulas  developed  by  Wheeler  [37]  give  2.54 
mm.  This  close  agreement  indicates  that  the  finite  element  method  is  correctly  predicting  the 
fields  at  the  interface  between  highly  contrasting  dielectrics,  even  in  the  presence  of  a  sharp 
conducting  edge. 
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7.3.2.  Microstrip  Meander  Line.  Since  the  modeling  method  accurately  accounts  for  all 
relevant  boundary  conditions  for  a  simple  microstrip  line,  it  should  be  obvious  that  it  can  correct¬ 
ly  predict  the  performance  of  most  passive  printed-circuit  RF  components.  A  microstrip  meander 
line  was  chosen  as  a  demonstration  case  since  measured  data  was  available.  The  circuit  dimen¬ 
sions  (in  /on)  and  the  finite  element  mesh  (substrate  cells  only)  are  shown  in  Figure  26. 

The  measurement  reference  planes  are  at  the  centers  of  the  pads  at  each  end  of  the  circuit, 
while  those  used  in  the  calculations  are  at  the  points  where  the  lines  begin  to  taper  from  75/zm 
down  to  25pm,  so  the  difference  must  be  accounted  for  when  comparing  the  data.  The  circuit 
layout  shows  the  location  of  via  holes  that  provide  a  ground  for  the  measurement  probes.  At  the 
contact  points,  the  75/im  center  conductor  and  50/xm-wide  slots  form  a  50Q  coplanar  waveguide 
(CPW).  Using  formulas  given  by  Rowe  and  Lao  [38],  the  effective  relative  permittivity  for  that 
transmission  line  structure  was  found  as  ceff=6.29.  The  measured  Soj  phase  was  corrected  by 
adding  360°  •  (75/*m)(ferf)l/Vx0.  The  calculated  S2i  phase  was  also  corrected  by  subtracting  the 
excess  phase  introduced  by  the  coaxial  connectors,  computed  from  a  straight  length  of  microstrip 
line  as  discussed  in  the  preceding  section.  Figure  27  compares  the  measurements  and  calculations 
over  the  frequency  range  from  1  to  26  GHz.  The  slight  discrepancy  in  the  phase  slope  is 
attributable  to  the  differing  reactances  of  CPW-microstrip  (measurement)  and  coax-microstrip 
(calculations)  transitions. 

7.4.  Importance  of  Higher  Order  Modes 

An  important  issue  is  whether  there  is  in  fact  any  benefit  to  using  higher-order  modes  in 
conjunction  with  a  finite  element  solution.  The  alternative  is  to  extend  the  finite  element  mesh 
into  the  connecting  waveguides  far  enough  that  any  higher  modes  excited  by  the  cavity’s  contents 
and/or  apertures  have  decayed  to  insignificance  |39[.  The  answer  depends  on  two  factors:  first, 
the  mode  excitation  strengths;  and  second,  their  attenuation  constants  in  the  waveguides,  which 
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(a)  (b) 


Figure  26.  Microstrip  Meander  Line:  (a)  Wafer  Metallization  Dimensions 
(in  pm);  and  (b)  Finite  Element  Mesh  for  Substrate 
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Figure  27.  Comparison  of  Measured  and  Calculated  S2 1  Phase  for 
Microstrip  Meander  Line 
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depends  mainly  on  the  ratio  of  the  frequency  to  the  mode  cutoff  frequency. 

Figure  28  shows  how  the  magnitude  of  Sjj  converges  with  the  number  of  modes  in  the 
one-quadrant  coaxial  iris  problem  (Figures  21  and  22).  The  most  important  modes  are  the  TEM 
and  the  TE(1  modes.  The  latter  is  excited  at  about  -10  dB,  and  its  cutoff  frequency  in  APC-7 
coax  is  about  20  GHz.  At  18  GHz  it  decays  by  29  dB  per  wavelength,  so  in  order  for  it  to  decay 
below  -30  dB  (regarded  as  negligible  for  purposes  of  S  parameter  calculation),  the  mesh  would 
need  to  extend  roughly  %X  into  the  waveguide  in  each  direction.  The  finite  element  mesh  and 
the  interaction  matrix  would  then  contain  an  additional  500-1500  edges  and  unknowns, 
respectively. 

On  the  other  hand,  the  two  microstrip  problems  do  not  excite  any  higher  order  modes 
above  the  -30  dB  level,  so  those  devices  could  be  accurately  modeled  using  only  the  lowest-order 
mode  and  still  terminating  the  mesh  at  the  waveguide  aperture.  Thus,  the  use  of  the  mode  sum 
as  a  continuity  condition  can  make  the  finite  element  solution  more  efficient,  but  the  improvement 


Figure  28.  Reflection  Coefficient  for  One-Quadrant  Coaxial  Waveguide  Iris 
Calculated  with  Varying  Numbers  of  Higher-Order  Modes 


57 


over  the  dominant-mode-only  approach  is  problem-dependent. 

7.5.  Summary 

In  addition  to  validating  the  essential  properties  of  the  electromagnetic  solution  approach, 
the  TWOPORT  code  provided  valuable  data  on  practicality:  typical  mesh  sizes  and  execution 
times.  The  data,  summarized  in  Table  I,  includes  all  of  the  test  cases  discussed  in  this  chapter, 
plus  three  others:  The  rectangular  and  circular  waveguide  cases  were  short  sections  of  waveguide 
used  to  verify  dominant  mode  propagation;  and  the  coax-to-rectangular  waveguide  transition 
consists  of  a  right-angle  metal  launcher  in  a  section  of  rectangular  waveguide  approximately  A/3 
long.  The  meshes  for  these  three  problems  are  shown  in  Chapter  VIII,  where  they  are  used  as 
cavity  array  elements. 


Table  I.  Mesh  Sizes,  Iterations  and  Matrix  Solve  Time  for  TWOPORT  Test  Cases 
(Solve  time  is  CPU  time  on  VAX®  8650  minicomputer) 


Test  Case 

Volume/A3 

Number 
Mesh  Cells 

Cells/A3 

Number  of 
Iterations 

Solve 

Time 

(min.) 

Rectangular  Waveguide 

.149 

348 

2067 

.9N 

0.7 

Circular  Waveguide 

.120 

342 

2850 

1.3N 

1.4 

Iris  in  Coaxial  Waveguide 

.0135 

847 

63X103 

.9N 

12.5 

Circular  Waveguide 

Mode  Converter 

.376 

1504 

4000 

0.5N 

43.5 

Coax-to-Rectangular 
Waveguide  Transition 

.294 

2442 

8294 

4.7N 

226. 

Microstrip  Meander  Line 

9.4xI05 

4190 

4.5xl07 

8.4N 

524. 

Microstrip  Transmission  Line 

5.7xl0'4 

4662 

8.2xl0b 

7. IN 

683. 

The  first  column  in  Table  I  gives  the  total  volume  in  cubic  wavelengths  for  each  problem. 
For  the  simpler  problems,  the  number  of  cells  is  determined  primarily  by  the  sampling 
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requirement  of  10  nodes  per  wavelength.  More  complicated  structures,  however,  require  more 
cells  in  order  to  capture  fine  details  of  geometry.  Column  4  gives  the  number  of  iterations  in 
terms  of  N,  the  number  of  mesh  edges,  that  were  required  to  obtain  convergence  of  the  residua! 
norm  to  .001  of  its  initial  value  (using  the  zero  vector  as  an  inital  guess).  In  all  cases,  the 
maximum  number  of  iterations  was  less  than  10  times  the  number  of  edges  (unknowns).  The 
matrix  solve  time  per  frequency  sample  is  given  in  the  last  column.  These  results  indicate  that 
practical  design  problems  can  be  accomplished  on  typical  minicomputers  and  workstations. 

This  chapter  has  validated  two  key  elements  of  the  solution  approach.  First,  it  showed 
that  FEM  can  correctly  account  for  boundary  conditions  typical  of  antenna  problems:  conductors, 
conductor  edges  and  dielectric  interfaces.  Second,  it  showed  the  validity  and  usefulness  of  the 
waveguide  mode  integral  equation  as  a  continuity  condition.  These  results  validated  not  only  the 
generic  concepts,  but  the  computer  code  implementation  as  well.  The  latter  was  especially 
important  since  the  TWOPORT  code  included  most  of  the  structure  and  modules  needed  for  the 
subsequent  cavity  array  problem  solution  discussed  in  the  next  chapter. 
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VIII.  Validation  -  Cavity  Array  Problem 


The  second  development  stage  replaced  the  integral  equation  for  the  outlet  waveguide 
with  the  periodic  integral  equation.  This  allows  validation  of  the  periodic  radiation  condition 
without  the  complexity  of  side-wall  periodicity  conditions.  It  also  provides  a  tool  useful  for 
analyzing  the  properties  of  radiators  such  as  open-ended  waveguides,  cavity-backed  slots  and 
multimode  horns.  This  chapter  describes  its  implementation  in  a  computer  code  and  presents 
comparisons  of  its  calculations  with  published  results.  Also  presented  are  two  specific  examples 
of  radiators  that  are  beyond  the  capability  of  previous  methods:  a  pyramidal  horn;  and  a  rectan¬ 
gular  waveguide  with  a  coaxial  transition  in  close  proximity  to  the  aperture. 

8.1.  Computer  Code  Implementation 

A  FORTRAN  program  named  CAVIARR  (cavity  array)  implements  the  solution  to  this 
generic  problem,  depicted  in  Figure  7.  Figure  29  is  an  outline  of  its  actions.  Note  that  all  of 
these  actions  up  to  step  IV. C.  are  essentially  the  same  as  in  program  TWOPORT.  The  interior 
finite  element  matrix  and  the  inlet  waveguide  submatrix  calculations  are  exactly  the  same,  since 
they  do  not  depend  on  the  scan  angle.  For  each  separate  scan  angle,  the  submatrix  due  to  the 
radiating  aperture  must  be  recalculated,  then  the  system  must  be  solved  again  for  the  field  vector 
E.  The  active  reflection  coefficient  calculation  is  the  same  as  the  calculation  for  Sn  in  TWO- 
PORT.  However,  instead  of  the  S2i  calculation,  CAVIARR  calculates  the  element  far  field  and 
transmission  coefficients  into  6  and  0  polarizations. 

8.2.  Waveguide  Arrays 

8.2.1.  Rectangular  Array.  The  scan-dependent  reflection  coefficients  for  several  open- 
ended  waveguide  arrays  are  available  from  publications  by  other  authors.  For  example,  the 
properties  of  the  rectangular  waveguide  array  shown  in  Figure  30a  were  first  presented  by  Dia- 
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I.  READ  INSTRUCTIONS  AND  OPTIONS 

II.  READ  PROBLEM  GEOMETRY 

III.  CREATE  EDGE-BASED  GEOMETRY 

IV.  FOR  EACH  FREQUENCY: 

A.  COMPUTE  TERMS  OF  SEE  ACCORDING  TO  (20),  (221 

B.  FOR  WAVEGUIDE  A: 

1.  COMPUTE  INCIDENT  CURRENT  VECTOR  IN  (32) 

2.  COMPUTE  TERMS  OF  S^  FROM  (34),  (35) 

C.  FOR  EACH  ANGLE: 

1.  COMPUTE  TERMS  OF  S^  FROM  (51) 

2.  SOLVE  (SEE  +  S0  )  E  =  Einc  FOR  E 

3.  COMPUTE  REFLECTION  COEFFICIENT  AND  MODE 
EXCITATION  COEFFICIENTS  FROM  (36) 

4.  COMPUTE  ELEMENT  FAR  FIELD  AND  TRANSMISSION 
COEFFICIENTS  FROM  (52)-(55) 

Figure  29.  Solution  Procedure  in  Program  CAVIARR 

mond  [40]  and  independently  confirmed  by  several  others  [29], [41], 

Figure  30b  is  the  finite  element  model  used  by  program  CAVIARR  -  simply  a  short 
section  of  waveguide.  The  nodes  on  the  shaded  surfaces  were  tagged  as  conductors,  while  those 
on  the  front  and  back  faces  were  tagged  as  radiating  aperture  and  waveguide  boundaries,  respec¬ 
tively. 

Figure  31  compares  Diamond’s  results  with  CAVIARR  computations.  The  interesting 
feature  of  this  test  case  is  the  scan  blindness  near  33°  in  the  H-plane  (the  4>  =  0  plane).  When  the 
array  scans  to  that  angle,  nearly  all  of  the  incident  power  is  reflected  back  into  the  waveguide. 
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The  CAVIARR  (HFEM)  calculations  used  waveguide  modes  up  to  m,n=2,3  and  Floquet  modes 
up  to  |m|, | n |  =5,5.  This  number  is  consistent  with  the  estimate  given  earlier  in  Section  5.3 
since  the  mesh  spacing  in  the  aperture  is  approximately  Xq/10.  In  fact,  |  m  | ,  |  n  |  <4  was  adequate 
to  ensure  convergence  of  the  active  reflection  coefficient  to  within  .  1  %  in  magnitude  and  1  °  in 
phase. 


8.2.2.  Rectangular  Array  with  Conducting  Iris.  A  straightforward  method  for  eliminat¬ 
ing  the  H-plane  scan  blindness  in  the  rectangular  waveguide  array  is  to  place  conducting  irises 
in  the  apertures.  For  example,  Lee  &  Jones  [41]  performed  a  multimode1  analysis  for  the  same 
array  lattice  and  waveguide  size  as  in  Figure  30,  except  that  a  thin  conductor  blocked  the  left  and 
right  14  of  each  aperture.  The  same  finite  element  mesh  as  before  (Figure  30b)  was  used  for  this 
problem,  except  that  the  nodes  associated  with  the  iris  were  tagged  as  separate  conductors. 
Figure  32  compares  the  CAVIARR  calculations  with  the  H-plane  element  gain  pattern  from  [41], 
The  scan  blindness  was  indeed  suppressed  by  the  addition  of  the  irises,  but  at  the  expense  of  a 
reduction  in  the  broadside  gain. 

The  close  agreement  between  the  CAVIARR  calculations  and  published  results  for  these 
two  test  cases  demonstrates  that  the  hybrid  periodic  integral  equation/finite  element  formulation 
is  valid.  It  also  demonstrates  that  its  implementation  in  the  computer  code  is  accurate. 

8.2.3.  Circular  Waveguide  Arrays.  The  lattice  geometry  for  an  array  of  circular  wave¬ 
guide  radiators  is  shown  in  Figure  33a.  The  incident  waveguide  field  is  in  the  dominant  TEn 
mode  and  polarized  parallel  to  the  y  axis.  Multimode  calculations  for  this  geometry  are  available 
from  Amitay  et.  al.  [29]. 

1  The  multimode  method  is  sometimes  referred  to  as  a  moment  method.  It  equates  the  transverse  fields 
on  the  two  sides  of  the  aperture  as  sums  over  modes  appropriate  to  each  region.  The  modes  of  one  waveguide 
are  used  as  testing  functions  in  order  to  generate  a  matrix  equation.  The  generic  technique  developed  for  wave¬ 
guide  discontinuities  was  adapted  to  the  phased  array  problem  using  Floquet  modes. 


63 


0  10  20  30  40  50  60 

SCAN  ANGLE  0  (deg.) 


Figure  34.  Circular  Waveguide  Array  Active  Reflection  Coefficient  vs.  Scan  Angle 

The  finite  element  mesh  used  for  CAVIARR  is  shown  in  Figure  33b.  Its  radius  was 
adjusted  so  that  the  sum  of  the  tetrahedron  volumes  was  the  same  as  a  perfect  circular  cylinder 
with  .343Aq  radius.  Figure  34  shows  the  finite  element  calculations  along  with  the  published 
multimode  results  for  both  E-  and  H-pIane  scanning  [29:276,280].  Again,  the  CAVIARR  code 
is  accurately  predicting  the  active  reflection  coefficient  at  all  scan  angles. 

A  method  for  suppressing  the  scan  blindness  in  circular  waveguide  arrays  is  to  use  dielec¬ 
tric-loaded  waveguides  that  can  be  packed  closer  together  since  their  radius  is  smaller  for  a  given 
wavelength.  In  the  following  test  case,  the  waveguide  radius  was  adjusted  so  that  eT'A  a  = 
.343Xo.  The  lattice  spacing  was  adjusted  so  that  dx/a  and  dy/a  are  constant.  The  finite  element 
model  was  simply  a  scaled  version  of  the  mesh  in  Figure  33b.  Figure  35  is  a  comparison  of  the 
loaded  and  unloaded  cases,  the  latter  with  er=4.1,  for  H-plane  scanning  [29:290], 

A  third  and  final  test  case  involving  circular  waveguide  arrays  also  includes  a  dielectric, 
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Figure  35.  Active  Reflection  Coefficient  vs.  Scan  Angle  -  Circular  Waveguide 
Array  with  and  without  Dielectric  Loading 

but  instead  of  completely  filling  the  waveguide,  it  is  only  a  short  plug  in  the  aperture  end,  flush 
with  the  ground  plane.  Amitay  et.  al.  give  results  for  E-plane  scanning  using  the  lattice  in  Figure 
33a,  and  .429Xo  length  plugs  with  er=  1.8  [29:293].  The  finite  element  model  for  this  problem 
was  a  longer  version  of  that  depicted  in  Figure  33b  -  it  was  6  mesh  cells  in  length,  in  order  to 
meet  the  A/10  sampling  requirement.  Figure  36  compares  those  results  for  the  no-dielectric  case. 
The  results  of  the  CAVIARR  and  multimode  computations  are  substantially  identical.  (There  are 
sources  of  error  in  both  approaches,  such  as  the  number  of  modes  and  the  number  of  integration 
points.  The  important  fact  is  that  both  methods  are  accurate  enough  to  allow  a  judgement  of 
whether  or  not  an  antenna  design  is  acceptable.)  This  demonstrates  the  important  capability  of 
the  hybrid  finite  element  method  for  modeling  arrays  that  have  dielectric  loading.  Unlike  the 
multimode  method,  it  is  not  restricted  to  homogeneous  dielectrics. 
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Figure  36.  Active  Reflection  Coefficient  vs.  Scan  Angie  -  Circular  Waveguide 
Array  with  and  without  Dielectric  Plugs 

8.3.  Pyramidal  Horn  Array 

The  multimode  method  was  the  only  rigorous  technique  available  for  addressing  cavity 
array  problems.  This  confined  the  available  solutions  to  structures  for  which  mode  sets  can  be 
constructed,  mostly  cylindrical  waveguides.  There  is,  however,  one  example  by  Amitay  &  Gans 
[42]  that  uses  a  variation  of  the  multimode  method  to  approximate  the  scanning  characteristics 
of  an  array  of  pyramidal  horns.  Their  technique  models  the  pyramidal  horn,  shown  in  Figure 
37a,  as  a  rectangular  waveguide  whose  dimensions  are  the  same  as  the  horn  mouth,  containing 
planes  at  various  z  locations  that  are  transparent  to  those  waveguide  modes  that  may  propagate, 
and  shorting  those  that  are  cut  off.  In  spite  of  the  approximate  nature  of  their  prediction  tech¬ 
nique,  they  obtained  fairly  good  agreement  with  measured  data. 

This  problem  is  a  fairly  stressing  case  for  the  hybrid  finite  element  method.  First  of  all, 
the  mesh,  depicted  in  Figure  37b,  has  approximately  6000  edges.  Furthermore,  the  large  unit 
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a  =  .905A, 


b  -  .452A 


dy-  1.19A. 

(a)  (b) 

Figure  37.  Pyramidal  Horn  Radiator:  (a)  Dimensions;  (b)  Tetrahedron  Mesh 

cell  size  (the  same  as  the  horn  mouth)  requires  an  unusually  large  number  of  Floquet  modes  (60 
in  kx  and  15  in  ky). 

Figure  38  is  a  comparison  of  the  CAVIARR  calculations  for  active  element  gain  with 
measured  data  from  [42],  The  scan  plane  is  0  =  90°,  for  which  the  co-polarized  field  is  the 
component.  The  gain  in  decibels  is  20  log10Ee.  The  interesting  feature  of  this  test  case  is  the 
scan  blindness  near  40°,  due  to  excitation  of  a  higher-order  waveguide  mode  at  the  aperture.  In 
contrast  to  the  waveguide  cases  discussed  earlier,  the  active  reflection  coefficient  does  not  become 
large  near  the  blindness  angle.  Since  there  is  more  than  one  propagating  Floquet  mode  in  this 
instance,  incident  power  that  is  not  transmitted  in  the  desired  direction  is  transmitted  into  another 
mode,  i.e.  a  grating  lobe. 
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8.4.  Coaxial-to-Rectangular  Waveguide  Launcher 


An  end-wall  transition,  or  "launcher,"  from  coaxial  to  rectangular  waveguide  is  shown 
in  Figure  39.  This  type  of  transition  hac  an  advantage  in  a  rectangular  waveguide  phased  array 
because  the  waveguides  are  packed  too  closely  to  use  a  broad-wall  launcher.  This  design,  due 
to  Tang  and  Wong  [43],  is  an  extension  of  the  coax  center  conductor  with  a  shorting  post  joining 
it  to  the  broad  wall.  The  probe  is  centered  in  height  and  offset  1/1  Oth  the  waveguide  width.  Its 
length  should  be  approximately  one  fourth  the  guide  wavelength  at  the  design  frequency.  For 
a  length  of  9mm,  the  frequency  response  is  shown  in  Figure  40  (computed  by  program  TWO- 
PORT). 

The  experimental  array  reported  in  [43]  had  a  long  waveguide  section  between  the  launch¬ 
ers  and  the  open  apertures,  and  higher-order  waveguide  modes  excited  by  either  one  would  not 
affect  the  other.  The  antenna’s  weight  and  size  would  both  benefit  if  the  launchers  were  placed 
as  close  as  possible  to  the  apertures,  but  then  their  mutual  interactions  are  not  negligible. 
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Figure  38.  Active  Element  Gain  vs.  E-Plane  Scan  Angle  for  Pyramidal 
Horn  Array  (Measured  Data  from  Amitay  &  Gans  1 42 ] ) 


REFLECTION  COEFFICIENT  MAGNITUDE 


(a)  (b) 

Figure  39.  End-Wall  Transition  from  Coax  to  Rectangular  Waveguide  (After  Tang  & 
Wong  [43]):  (a)  Geometry;  (b)  Cutaway  of  Tetrahedron  Mesh 
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Figure  40.  Reflection  Coefficient  (S  j  j )  for  Coax-Rectangular  Launcher. 
Probe  Length  =  9mm;  Post  Height  =  5mm;  X-BanJ  Waveguide 
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Considering  again  the  rectangular  waveguide  test  case  discussed  in  Section  8.2.1,  suppose  the 
waveguide  elements  are  only  10mm  deep,  with  a  launcher  extending  to  within  2.5mm  of  the 
aperture.  The  predicted  active  reflection  coefficient  vs.  angle  is  shown  in  Figure  41,  compared 
to  the  original  case,  where  the  array  is  fed  by  semi-infinite  rectangular  waveguides  (both  calcula¬ 
tions  by  CAVIARR).  Evidently,  the  launcher  in  close  proximity  to  the  waveguide  opening 
prevents  the  formation  of  the  higher-order  waveguide  mode  that  is  responsible  for  the  scan 
blindness  condition,  which  was  the  purpose  of  the  conducting  irises  discussed  by  Lee  &  Jones. 
Furthermore,  the  interactions  between  the  launcher  and  the  waveguide  aperture  do  not  generate 
any  additional  resonance  effects,  so  there  is  no  need  to  separate  the  two  by  a  long  length  of 
waveguide.  Thus,  the  array  element  consisting  of  a  coaxial  launcher  in  a  short  waveguide  section 
is  a  better  solution  than  was  previously  available  because  it  achieves  the  same  result  with  a 
smaller  and  simpler  structure.  Although  rectangular  waveguide  arrays  are  outdated,  this  is 
nonetheless  an  illustration  that:  (a)  the  hybrid  finite  element  method  can  be  used  to  solve  practical 


Figure  41.  Active  Reflection  Coefficient  for  Rectangular  Waveguide  Array  with 
Coaxial  Launcher  and  Rectangular  Waveguide  Feeds 
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problems  for  which  accurate  methods  were  not  available  before;  and  (b)  the  cavity  array  solution 
has  important  applications  by  itself  even  though  it,  like  TWOPORT,  was  an  intermediate  step  in 
the  code  development.  The  following  chapter  describes  the  implementation  and  results  from  the 
final  step,  which  implements  the  periodicity  conditions  at  unit  cell  side  walls  for  general  array 
radiators. 
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IX.  Validation  -  General  Array  Problem 


The  third  and  final  stage  of  computer  code  development  and  validation  implemented  the 
periodicity  condition  for  unit  cell  side  walls,  removing  the  requirement  that  the  radiators  be 
separated  by  conducting  walls.  Initial  validations  were  performed  for  open  waveguide  arrays 
using  the  previous  chapter’s  results.  Test  cases  for  arrays  of  microstrip  patches  and  monopole 
arrays  were  used  for  further  validation,  the  latter  including  a  hardware  experiment.  Last,  perfor¬ 
mance  predictions  for  flared  notch  and  printed  dipole  arrays  are  presented  to  demonstrate  the 
method’s  versatility. 

9.1.  Computer  Code  Implementation 

A  FORTRAN  program  named  PARANA  (phased  array  antenna  analysis)  implements  the 
general  array  problem  solution.  Figure  42  is  an  outline  of  the  program’s  actions.  Up  to  step 
IV.C.  it  is  essentially  identical  to  CAVIARR,  except  that  a  listing  of  image  edges  is  created  in 
step  III.  These  identify  the  -x  and  -y  boundary  counterparts  for  each  edge  on  the  +x  and  +y 
unit  cell  walls.  This  places  an  additional  requirement  on  the  geometry  file  -  the  unit  cell  bound¬ 
ary  nodes  must  be  tagged  so  that  these  edges  can  be  identified.  In  addition,  it  requires  the  mesh 
generator  to  create  the  grid  in  such  a  way  that  the  surface  mesh  on  opposing  faces  is  identical. 
Although  the  CAD  software  (I-DEAS™  [44])  does  not  have  any  provision  for  enforcing  this 
requirement,  it  results  naturally  when  the  mesh  areas  and  constituent  curves  are  defined  in 
consistent  and  logical  order. 

In  step  IV.C.  1 .,  the  calculation  for  the  radiation  boundary  terms  of  is  slightly  modi¬ 
fied  to  include  special  handling  for  -x  and  -y  boundary  edges,  as  discussed  in  Section  6.4.  Step 
IV.C. 2.  is  a  straightforward  implementation  of  the  algorithm  shown  in  Figure  17.  One  of  its 
effects  is  to  zero  out  the  rows  and  columns  of  SEE  corresponding  to  +x  and  +y  boundary  edges. 
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I.  READ  INSTRUCTIONS  AND  OPTIONS 

II.  READ  PROBLEM  GEOMETRY 

III.  CREATE  EDGE-BASED  GEOMETRY  AND  LIST  OF  IMAGE  EDGES 

IV.  FOR  EACH  FREQUENCY: 

A.  COMPUTE  TERMS  OF  SEE  ACCORDING  TO  (20),  (22) 

B.  FOR  WAVEGUIDE  A: 

1.  COMPUTE  INCIDENT  CURRENT  VECTOR  FROM  (33) 

2.  COMPUTE  TERMS  OF  S^  FROM  (34) 

C.  FOR  EACH  ANGLE: 

1.  COMPUTE  TERMS  OF  S^  FROM  (51)  USING  OVERLAP  ELEMENTS 

2.  IMPOSE  PERIODICITY  CONDITION  ON  SEE  ACCORDING  TO  FIG.  17 

3.  ELIMINATE  ZERO  ROWS  &  COLUMNS  FROM  (SEE  +  S^  ) 

4.  SOLVE  (SEE  +  S^  )  E  =  Einc  FOR  E 

5.  RESTORE  BOUNDARY  EDGES 

6.  COMPUTE  REFLECTION  COEFFICIENT  AND  MODE 
EXCITATION  COEFFICIENTS  FROM  (36) 

7.  COMPUTE  ELEMENT  FAR  FIELD  AND  TRANSMISSION 
COEFFICIENTS  FROM  (52)-(55) 

Figure  42.  Solution  Procedure  in  Program  PARANA 

Step  IV. C.  I.  had  a  similar  effect  on  S^.  Before  starting  the  matrix  solution,  the  matrix  is 
compressed  to  eliminate  zero  rows  and  columns.  After  the  matrix  is  solved,  the  periodicity 
conditions  are  used  to  solve  for  the  field  along  the  +x  and  +y  boundary  edges. 

9.2.  Wa  ve guide  A  rrays 

The  PARANA  code  is  also  capable  of  modeling  open-ended  waveguide  arrays,  although 
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less  efficiently  than  CAVIARR.  A  comparison  between  the  two  was  a  means  for  verifying  the 
periodic  boundary  condition  while  keeping  all  other  details  of  the  problem  and  solution  the  same. 
Figure  43  is  the  tetrahedron  mesh  used  as  the  input  to  PARANA.  It  is  a  shallow  slice  of  the  unit 
cell  in  free  space  outside  the  waveguide,  above  the  ground  plane.  The  shaded  area  identifies 
those  nodes  and  cells  bordering  on  the  ground  plane.  Contrast  this  to  the  CAVIARR  geometry 
model  (Figure  33b),  which  is  a  slice  of  waveguide  below  the  ground  plane. 

The  test  case  parameters  were:  a  =  JXq  ,  dx  =  .88Xo  ;  dy  =  .733XQ  and  y  =  60° 
(equilateral  lattice).  Figure  44  compares  the  computed  results  from  the  two  codes  for  scanning 
in  the  H-plane  (x-z  plane).  The  fact  that  they  are  essentially  identical  is  evidence  that  the  periodic 
boundary  condition  is  working  as  proposed.  The  results  for  E-plane  scanning  were  similar,  with 
the  same  degree  of  agreement  between  the  two  codes. 

A  second  waveguide  test  case  was  a  rectangular  lattice  of  rectangular  waveguides.  The 
lattice  dimensions  were  the  same  as  the  waveguide  size,  i.e.:  a  =  dx  =  23mm  and  b  =  dy  = 
10mm.  The  finite  element  mesh  was  similar  to  Figure  30  except  that  there  were  three  layers  of 
cells,  each  with  the  same  thickness.  The  perimeter  edges  of  the  first  layer  were  tagged  as  con- 


Figure  43.  Finite  Element  Mesh  for  a  Skewed-Lattice,  Circular  Waveguide  Array  Unit  Cell 
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Figure  44.  Circular  Waveguide  Array  Active  Reflection  Coefficient  -  Comparison  of  Results 
Using  Cavity  Array  (CAVIARR)  and  General  Array  (PARANA)  Models 

ductors  to  indicate  that  the  mesh  extends  part  way  into  the  waveguide.  The  perimeter  edges  of 
the  other  two  layers  were  tagged  as  unit  cell  walls.  The  calculations  by  PARANA  for  10  GHz, 
shown  in  Figure  45,  are  again  essentially  the  same  as  the  CAVIARR  results.  This  test  case 
illustrates  some  potentially  important  flexibilities  of  the  implementation:  the  mesh  may  extend 
into  the  feed  waveguide;  and  the  feed  waveguide  may  adjoin  the  unit  cell  boundary. 

9.3.  Microstrip  Patch  Array 

The  first  of  four  demonstration  cases  uses  a  simple  rectangular  microstrip  patch  on  a 
substrate  that  is  thick  enough  to  support  a  surface  wave  in  the  dielectric.  The  geometry,  shown 
in  Figure  46,  has  the  dimensions  (in  wavelengths)  used  by  Pozar  in  MoM  calculations  [45].  The 
substrate  has  relative  permittivity  of  12.8  to  represent  GaAs.  Since  the  input  impedance  of  the 
patch  is  very  low  (less  than  50)  due  to  the  substrate  thickness,  it  presents  a  large  mismatch  to  the 
coaxial  feed.  Consistent  with  [45],  the  active  reflection  coefficient  is  normalized  using: 
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ACTIVE  REFLECTION  COEFFICIENT 


SCAN  ANGLE  (deg) 

Figure  45.  Rectangular  Waveguide  Array  Active  Reflection  -  Comparison  of  Results 
Using  Cavity  Array  (CAVIARR)  and  General  Array  (PARANA)  Models 


Figure  46.  Microstrip  Patch  Radiator  (dimensions  in  wavelengths) 
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(65) 


JVM)  = 

zin(d,4>)+z*{0,0) 

This  measure  discounts  the  effect  of  the  feed  impedance  mismatch  to  more  clearly  reveal  the 
effects  of  scanning.  Figure  47  compares  the  PARANA  (HFEM)  and  moment  method  calcula¬ 
tions.  The  agreement  is  within  approximately  7  % ,  with  the  greatest  difference  at  the  angle  where 
there  is  a  scan  blindness  due  to  a  surface  wave  in  the  dielectric  slab.  This  discrepancy  may  be 
due  to  the  feed  modeling  ([45]  used  an  idealized  probe  feed).  The  fact  that  the  PARANA  code 
is  predicting  the  existence  of  the  trapped  surface  wave  is  a  further  confirmation  that  the  periodic 
boundary  conditions  are  effective  and  correctly  implemented. 

9.4.  Clad  Monopole  Array  Experiment 

A  candidate  antenna  design  for  a  space-based  radar,  proposed  by  Fenn,  was  a  planar 
array  of  monopoles  above  a  ground  plane  [46].  Each  monopole  is  simply  an  extension  of  a 


Figure  47.  Active  Reflection  Coefficient  for  Microstrip  Patch  Array.  E-Plane  Scan 
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coaxial  waveguide’s  center  conductor.  Fenn’s  MoM  analysis  was  adequate  for  the  bare  mono¬ 
pole,  but  not  for  one  with  a  dielectric  sheath,  or  cladding.  The  sheath,  which  can  be  simply  an 
extension  of  the  coax  insulator,  is  an  important  enhancement  to  the  monopole  design  because  it 
increases  the  bandwidth  substantially. 

9.4.1.  Initial  Validation.  For  a  first  validation  check,  PARANA  calculations  were 
compared  to  Fenn’s  calculations  for  .25\)-Iong  monopoles  arrayed  in  a  .5Xo  x  ,5Xq  lattice. 
Figure  48  shows  the  active  reflection  coefficient  vs.  scan  angle.  The  essential  scanning  charac¬ 
teristics  are  verified,  although  there  is  some  discrepancy  at  the  angles  where  the  reflection  coef¬ 
ficient  is  very  small.  This  is  at  least  partly  due  to  a  simplification  in  Fenn’s  MoM  model:  It 
used  an  assumed  form  of  current  distribution  (piecewise  sinusoidal)  on  the  monopole  that  is  not 
allowed  to  change  with  scan  angle.  Herper  and  Hessel  used  the  same  simplification,  and  showed 
a  large  disparity  between  calculated  and  measured  results  in  the  vicinity  of  60°  [47J. 

9.4.2.  Bandwidth  Enhancement  with  Cladding.  Figure  49  shows  one  half  of  the  tetrahe¬ 
dron  mesh  used  by  PARANA  for  the  following  test  case:  The  monopole  (15.875mm  long)  is 
represented  as  a  void  extending  through  the  mesh  from  top  to  bottom.  Shading  has  been  added 
to  the  figure  to  identify  the  cells  comprising  the  dielectric  cladding  (Teflon,  er=2.1).  The 
remaining  mesh  cells  are  free  space.  The  lattice  is  triangular  with  dx  =  36mm  and  dy  =  31mm 
(nearly  equilateral).  This  interelement  spacing  allows  scanning  over  the  0<9O°  hemisphere 
without  grating  lobes  for  frequencies  up  to  4.8  GHz. 

Figure  50  shows  contour  plots  of  |RJ  vs.  scan  ang'e  and  frequency  for  clad  and  unclad 
monopoles.  These  calculations  are  for  the  (j>  =  90°  scan  plane,  but  similar  results  would  be  ob¬ 
served  for  any  scan  plane  due  to  the  equilateral  lattice.  Both  arrays  achieve  an  acceptable  reflec¬ 
tion  coefficient  (|Ra|  ^.33,  or  VSWR^2:1)  over  only  a  very  limited  range  of  scan  angles  near 


79 


Figure  48.  Reflection  Coefficient  vs.  Scan  Angle  (</>=0  plane)  for  a  Monopole  Array 


Figure  49.  Finite  Element  Mesh  for  a  Clad  Monopole  in  a  Triangular  Lattice 


81 


I 


FREQUENCY  (GHz)  FREQUENCY  (GHz) 


6=60°.  The  clad  monopole  array  appears  to  have  nearly  twice  the  bandwidth  of  the  unclad  array 
at  6= 60°.  This  simple  example  case  demonstrates  the  usefulness  of  this  code  for  investigating 
radiator  designs  -  the  improved  (clad)  radiator  could  not  be  simulated  with  existing  method  of 
moments  codes. 

9.4.3.  Experiment.  A  121-element  clad  monopole  array  was  fabricated  using  the  same 
lattice  parameters  (dx  =  36mm,  dy=31mm)  in  an  isosceles  lattice.  The  monopole  element 
simply  an  unmodified  SMA  connector.  When  installed  in  the  l/8"-thick  Aluminum  backplane, 
the  conducting  monopole  and  dielectric  cladding  extend  14.7mm  and  11.8mm,  respectively, 
above  the  ground  plane. 

The  array  layout  is  shown  in  Figure  51 .  The  measurement  procedure  consists  of  connect- 


Figure  51.  Experimental  Array  Geometry  used  for  Coupling  Measurements 
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ing  the  center  element  to  port  #1  of  a  network  analyzer,  then  connecting  each  other  element  in 
turn  to  port  #2  and  measuring  the  coupling  coefficient  Cinn. 

The  active  reflection  coefficient  as  a  function  of  angle  is  given  by 


Ra(6,4>)  ~ 


cmn  e  e 


m  -  -oo  n  =  ~o° 


Once  a  finite  number  of  C^’s  are  found  by  measurement,  an  approximation  to  Ra  may  be 
computed  for  all  angles.  Figure  52  compares  measured  results  and  PARANA  calculations.  The 
oscillation  (with  angle)  of  the  measured  data  about  the  actual  Ra  was  expected  due  to  the  finite 
array  size. 

These  results,  combined  with  the  earlier  waveguide  and  microstrip  patch  results  provided 
the  validation  for  the  PARANA  code.  The  final  two  radiator  designs  discussed  next  are  attempts 
to  exploit  the  code  to  assess  the  properties  of  radiators  for  which  results  are  not  available  by  other 
computational  methods. 


9. 5.  Printed  Dipole  Radiator 

9.5.1.  Element  Design.  The  design  for  the  dipole  element  shown  earlier  in  Figure  2 
follows  general  guidelines  given  by  Edward  &  Rees  [6].  Its  height  is  . 25\ its  overall  length 
is  .5A o  and  its  arms  are  .OSXq  wide.  The  design  center  frequency  is  4.8  GHz,  giving  \0= 
62.5mm.  The  substrate  material  chosen  for  this  case  is  50  mils  thick  with  relative  permittivity 
of  10  because  such  material  is  readily  available  and  it  provides  a  good  scaled  representation  of 
100/un  GaAs  (Gallium  Arsenide)  at  millimeter  wave  frequencies. 

The  actual  dimensions  used  for  this  test  case  are  shown  in  Figure  53.  Figure  54  is  an 
exploded  view  of  the  tetrahedron  mesh  used  as  the  input  to  PARANA,  showing  the  three  mesh 
regions  corresponding  to  the  substrate  and  the  two  air  regions  on  each  side  filling  the  unit  cell. 
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Figure  52.  Measured  and  Computed  Active  Reflection  Coefficient  vs.  Angle 
for  Clad  Monopole  Array  Experiment 


The  array  lattice  is  square,  with  ,5X0  inter-element  spacing.  Note  that  the  mesh  is  denser  in  the 
dielectric  slab  (by  a  factor  of  er1/2  ),  and  gradually  relaxes  going  out  towards  the  sides  of  the  unit 
cell.  The  dipole  and  balun  are  metallized  (or  photoetched)  on  opposite  sides  of  the  substrate. 
The  slot  in  the  dipole  center  divides  that  part  of  the  structure  into  two  coupled  microstrip  lines. 
They  are  to  be  driven  180°  out  of  phase  by  the  balun.  The  dipole  wid*h  is  chosen  as  three  times 
the  balun  width  so  that  it  provides  a  ground  plane  for  the  microstrip  lines  comprising  the  balun. 
In  the  design  shown  in  Figure  53,  the  first  arm  of  the  balun  is  a  linear  taper  from  microstrip  line 
widths  corresponding  to  50Q  (the  same  as  the  coaxial  input)  to  8011  (the  narrow  end).  The  second 
and  third  balun  arms  must  have  the  same  characteristic  impedance  as  the  coupled  microstrips,  and 
values  below  80Q  require  extremely  narrow  slots  that  are  difficult  to  construct  given  the  toleranc¬ 
es  of  photolithographic  processes.  The  slot  length  is  approximately  1/4  guide  wavelength  from 
its  closed  end  to  where  it  crosses  underneath  the  microstrip  line.  Similarly,  the  microstrip  line 
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Figure  53.  Printed  Dipole  Radiator  Design 


Figure  54.  Exploded  Finite  Element  Mesh  for  Printed  Dipole  Radiator 


length  is  1/4  guide  wavelength  from  the  slot  to  its  open  end,  but  it  is  reduced  by  an  effective 
length  using  an  approximate  formula  due  to  Hammerstad  (48]  (see  also  [49:190]). 

9.5.2.  Calculations.  Initial  calculations  of  Ra  vs.  frequency  revealed  that  the  design  in 
Figure  53  was  a  poor  radiator,  with  jRJ  greater  than  .8  at  the  design  center  frequency  (4.8 
GHz).  Slightly  better  results  were  obtained  by  scaling  the  height  by  .8  so  that  the  overall  height 
was  12.5mm  (.20^)  and  then  reducing  the  dipole  width  from  .4Xq  to  .33Xq.  Figure  55  shows 
|  Ra  |  vs.  frequency  for  both  dipole  widths  with  the  12.5mm  height.  Unfortunately,  this  printed 
dipole  has  not  been  tested  as  an  isolated  radiator  on  high-permittivity  substrates,  so  it  is  not 
known  how  much  of  the  mismatch  is  due  to  array  effects  vs.  feed  line  mismatches.  Nonetheless, 
the  scanning  properties  may  still  be  evaluated,  normalizing  the  active  reflection  coefficient  using 
(65).  The  results  for  4.8  GHz,  shown  in  Figure  56,  indicate  that  there  are  no  scan  blindnesses 
in  either  the  E-  or  H-plane  (the  E  plane  is  the  plane  containing  the  substrate). 


Figure  55.  Computed  Active  Reflection  Coefficient  at  Broadside  Scan  for  Reduced-Height 
(12.5mm)  Printed  Dipole  for  Two  Dipole  Lengths 
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Figure  56.  Computed  Active  Reflection  Coefficient  (Normalized)  vs.  Scan  Angle 
for  Reduced  Height  (12.5mm)  Printed  Dipole,  4.8  GHz 


9.6.  Flared  Notch  Radiator. 

9.6.1.  Element  Design.  The  basic  idea  behind  the  flared  notch  radiator  shown  earlier 
in  Figure  1  is  a  slotline,  gradually  opening  out  to  provide  a  tapered  impedance  match  to  free 
space.  There  do  not  appear  to  be  any  specific  design  rules  [50],  but  generally,  the  longer  the 
flare,  the  greater  the  bandwidth.  For  purposes  of  this  study,  the  exponential  flare  given  by 
Choung  &  Chen  was  selected  [51]. 

Figure  1  showed  the  slotline  being  fed  by  a  microstrip  line,  which  is  in  turn  fed  by  a 
coaxial  cable.  This  is  a  form  of  balun  (the  same  arrangement  used  for  the  dipole  in  the  preceding 
section)  matching  the  balanced  coax  to  the  unbalanced  slotline.  An  alternative  balun  design  is 
based  on  a  new  coplanar  waveguide  (CPW)-to-slotline  transition  that  terminates  one  side  of  the 
CPW  in  a  broadband  open  circuit  [52],  This  design  has  the  advantage  that  only  one  side  of  the 
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substrate  is  metallized,  reducing  the  number  of  steps  in  fabrication  and  eliminating  the  possibly 
of  registration  error.  The  metallization  pattern  for  a  single  flared  notch  radiator  is  the  shaded 
area  in  Figure  57.  The  test  case  used  a  50-mil  thick  substrate  with  relative  permittivity  of  6.C. 
The  flare  length  and  mouth  width  are  33.3mm  and  30.0mm,  respectively.  The  flare  shape  is 
given  by 

w(z)  =  w0expjiln 

where  w0  and  wm  are  the  widths  at  the  slotline  an  the  mouth  and  L  is  the  flare  length. 


9.6.2.  Array  Performance  Calculations.  One  of  the  geometry  models  used  as  the  input 
to  PARANA  is  shown  in  Figure  58.  As  was  the  case  with  the  dipole  element,  the  mesh  is  denser 
in  the  substrate  than  in  the  free  space  regions.  The  unit  cell  size  is  dx=36mm  and  dy  =  34mm 
(rectangular  lattice).  A  second  model  (not  shown)  used  an  equilateral  triangular  lattice  with 
dx=62mm,  dy  =  36mm  and  7=60°.  Both  of  them  used  feed  waveguide  dimensions  corresponding 
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OPEN  CIRCUIT 


Figure  57.  Flared  Notch  Element  Design 
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GROUND 
PLANE  FACE 


SUBSTRATE  SURFACE  MESH 


(a) 


(b) 


Figure  58.  Finite  Element  Mesh  for  Flared  Notch  Radiator  (Rectangular  Lattice): 
(a)  Unit  Cell  Showing  Coax  Aperture;  (b)  Substrate  Surface  Mesh 


to  APC-3.5  coax  (a=  ,75mm,  b=  1.75mm). 

The  broadside  (0o=O)  active  reflection  coefficient  is  shown  as  a  function  of  frequency  in 
Figure  59.  These  indicate  that  the  radiator  is  capable  of  very  broad  bandwidth,  but  the  skewed- 
lattice  array  has  a  resonance  effect  that  causes  a  blindness  near  4.25  GHz.  The  rectangular-lattice 
array  is  a  very  promising  design,  since  its  predicted  bandwidth  of  greater  than  50%  is  difficult 
to  achieve  in  an  array. 

The  active  reflection  coefficient  vs.  scan  angle  for  the  rectangular-lattice  array  is  shown 
in  Figure  60.  From  the  low  (3GHz)  to  the  center  (4GHz)  of  the  frequency  range,  it  behaves 
well,  but  at  the  high  end  (5GHz)  it  displays  blindnesses  in  both  scan  planes  (the  E  plane  is  the 
plane  containing  the  substrate). 
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Figure  59.  Active  Reflection  Coefficient  vs.  Frequency  for  Printed  Flared  Notch  Arrays 

9. 7.  Summary 

The  results  of  this  chapter  have  proven  the  validity  of  the  essential  feature  that  makes  the 
hybrid  finite  element  method  applicable  to  infinite  array  analysis:  the  periodic  boundary  condi¬ 
tion  implemented  by  "wrapping"  opposing  unit  cell  edges  onto  each  other  with  an  appropriate 
phase  shift.  It  was  successful  for  both  rectangular  and  trapezoidal  unit  cells,  the  latter  applying 
to  skewed  array  lattices.  The  microstrip  patch  array  test  case  showed  that  it  correctly  predicts 
the  behavior  of  surfaces  waves  in  a  dielectric  layer  on  a  ground  plane.  The  monopole  array  test 
case  further  demonstrated  the  method’s  ability  to  deal  with  inhomogeneous  dielectrics.  Most 
importantly,  the  same  computer  program  with  no  changes  whatsoever  executed  the  computations 
for  every  one  of  the  seven  separate  array/radiator  geometries  discussed  in  this  chapter.  The  only 
things  that  changed  were  the  finite  element  model  created  in  I-DEAS™  and  the  user  instructions 
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identifying  the  feed  waveguide  type  and  location  and  the  array  lattice  parameters.  Thus  it  is 
demonstrated  that  the  finite  element  method  has  unparalleled  versatility  for  phased  array  analysis. 

The  question  of  efficiency  is  addressed  below  in  Table  II,  which  summarizes  the  time  and 
storage  requirements  for  most  of  the  test  cases  discussed  in  this  chapter.  The  total  run  time 
(tabulated  figures  are  for  a  VAX®  4400  mincomputer  whose  performance  is  rated  at  about  17 
MIPS)  is  dominated  by  the  matrix  solve  time,  which  was  as  high  as  7  hours  per  point  for  the 
printed  dipole  (the  most  time-consuming  case).  The  most  time-consuming  part  of  the  matrix  fill 
is  the  calculation  of  the  radiating  aperture  terms.  It  is  highest  for  the  microstrip  patch  because 
that  case  had  relatively  more  edges  in  the  aperture  than  did  the  other  cases.  The  fact  that  these 
computations  could  be  performed  on  a  typical  minicomputer  are  encouraging,  although  whether 
or  not  one  regards  them  as  "efficient"  depends  on  the  difficulty  of  solving  the  design  problem  by 
other  means. 


Table  II.  Mesh  Size,  Execution  Time  and  Matrix  Storage  for  PARANA  Test  Cases 


Test 

Case 

Unit  Cell 
Vol.  (X3) 

Mesh 

Cells 

Mesh 

Edges 

Iter¬ 

ations 

Fill 

Time 

(min.) 

Solve 

Time 

(min.) 

Matrix 

Size 

(Mbytes) 

Microstrip  Patch 

.015 

4291 

4862 

2500 

44.5 

198 

6.2 

Clad  Monopole 

.063 

3092 

3487 

19000 

7.1 

201 

0.7 

Flared  Notch  with 
Triangular  Lattice 

.17 

4443 

4778 

18000 

8.7 

124 

0.8 

Flared  Notch  with 
Rectangular  Lattice 

.19 

4845 

5271 

23500 

9.5 

196 

0.9 

Printed  Dipole 

.05 

6659 

7447 

39000 

20.6 

428 

1.6 
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X.  Conclusions  and  Recommendations 


This  research  project  began  with  a  novel  concept  for  modeling  infinite  phased  arrays  and 
concluded  with  a  demonstration  of  its  capability.  The  work  in  between  involved  the  entire 
process  of  electromagnetic  predictive  code  development:  casting  the  physical  problem  as  a 
mathematical  boundary  value  problem;  mapping  the  latter  to  a  linear  system  using  finite  element 
and  moment  methods;  designing,  writing  and  troubleshooting  the  general-purpose  computer 
program;  adapting  the  commercial  software  for  geometry  generation;  and  finally,  validating  the 
code.  Many  difficulties  needed  to  be  overcome,  both  expected  and  otherwise;  yet,  other  expected 
problems  proved  inconsequential.  This  chapter  attempts  to  summarize  those  findings  and  to 
assess  the  implications  of  the  results  to  future  work. 

10. 1.  Conclusions 

10.1.1.  Theory  and  Formulation.  The  successes  of  other  researchers  in  applying  the 
finite  element  method  to  time-harmonic  electromagnetic  problems  was  reason  to  believe  that  it 
would  also  succeed  for  the  phased  array  problem.  Nonetheless  there  are  always  doubts  surround¬ 
ing  any  such  implementation  given  that  the  problem  does  not  have  the  properties  of  self-adjoint¬ 
ness  or  positive  definiteness.  The  attendant  risk  that  the  iterative  matrix  solver  might  converge 
to  a  false  solution,  or  not  converge  at  all,  did  not  materialize.  The  conclusion  is  that  the  weak 
form  and  Galerkin’s  method  are  appropriate  to  this  class  of  problems. 

There  have  not  been  any  cases  that  would  indicate  instability  or  non-uniqueness,  which 
are  always  concerns  when  the  properties  of  fields  in  two  or  more  regions  must  be  met.  A  unique 
solution  generally  requires  that  continuity  of  tangential  electric  and  magnetic  fields  must  be 
independently  enforced.  The  fact  that  this  problem  involves  boundaries  that  are  planar  and 
exterior  solutions  in  terms  of  discrete  modes  whose  tangential  field  components  are  dependent 
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evidently  circumvents  that  requirement. 

10.1.2.  Implementation.  Three-dimensional  finite  element  problems  involve  very  large 
systems  of  equatio.  s.  As  seen  repeatedly  in  the  previous  four  chapters,  even  relatively  small 
devices  result  in  grids  with  thousands  of  edges,  usually  because  of  the  need  to  capture  fine  d  'tails 
in  the  geometry.  The  uncertainty  of  whether  the  available  computers  could  solve  the  resulting 
matrices  in  a  reasonable  time  has  been  resolved  and  the  validation  tests  indicate  that  most  practi¬ 
cal  array  radiators  may  be  analyzed  using  typicJ  minicomputers. 

The  three  codes  TWOPORT,  CAVIARk  and  PARANA  are  geometry-independent,  within 
the  constraints  of  the  generic  classes  of  problems  they  are  designed  to  solve.  They  exploit  the 
commercial  CAD  software  that  benefits  from  decades  of  research  oriented  towards  mechanical 
engineering  applications.  Two  observations  regarding  that  software  are:  (a)  it  is  capable  of 
generating  grids  suitable  for  phased  array  analysis;  and  (b),  of  much  greater  significance,  it 
removes  the  geometry  generation  from  the  electromagnetic  problem  so  that  the  codes  may  encom¬ 
pass  far  broader  problem  classes  than  previously  attem  pted. 

10. 1.3.  Validation.  The  extensive  set  of  validation  cases  that  were  used  to  test  the  three 
computer  codes  demonstrate  the  effectiveness  of  the  key  elements  of  the  solution  approach.  The 
interior  finite  element  solution  accurately  incorporates  electromagnetic  boundary  conditions  at 
perfect  conductors  (including  sharp  edges)  and  dielectric  interfaces.  It  is  evidently  free  of  spuri¬ 
ous,  non-physical  solutions,  indicating  that  the  divergence  condition  is  also  satisfied. 

The  waveguide  mode  integral  equation  was  shown  to  be  an  effective  mechanism  for 
enforcing  field  continuity  at  waveguide  apertures.  Its  payoff  is  in  providing  an  accurate  means 
for  modeling  antenna  feed  structures. 

The  hybridization  of  the  periodic  integral  equation  with  the  finite  element  solution  was 
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shown  to  correctly  enforce  the  radiation  condition  above  the  infinite  array.  Its  implementation 
ion  in  the  CAVIARR  code  allowed  the  correct  prediction  of  active  reflection  coefficient  for  a 
variety  of  radiator  types,  including  cases  involving  scan  blindnesses. 

The  implementation  of  a  periodic  boundary  condition  on  a  three-dimensional  finite 
element  problem  is  evidently  a  first.  The  success  of  the  final  stage  of  this  work  hinged  entirely 
on  that  single  unproven  algorithm  and  the  question  of  whether  or  not  it  would  conflict  with  the 
periodic  integral  equation.  Thus,  the  most  important  finding  is  that  the  algorithm  involving 
boundary  "wrapping"  works  as  proposed,  and  complements  the  periodic  radiation  condition. 

10.2.  Recommendations 

The  recommendations  fall  into  three  broad  categories.  The  first  deals  with  hardware 
experiments  needed  to  demonstrate  the  potential  for  design ,  in  contrast  tt.  the  analysis  that  was 
the  main  subject  of  this  project.  The  initial  designs  for  the  printed  dipole  and  flared  notch 
radiators  were  a  first  cut,  and  their  performance  leaves  much  to  be  desired.  Further  experiments 
will  be  necessary  to:  first,  perfect  the  design  for  single,  isolated  radiators;  second,  identify  the 
geometry  parameters  that  influence  their  impedance  match  in  the  array  environment;  and  finally, 
use  the  codes  to  optimize  those  parameters  for  best  performance  over  some  specified  range  of 
scan  angles  and  frequency. 

The  second  category  deals  with  possible  improvements  and  enhancements  to  the  computer 
codes.  For  example,  the  TWOPORT  code  could  easily  be  extended  to  deal  with  multi-port  RF 
devices.  The  inclusion  of  that  feature  in  PARANA  would  allow  the  simulation  of  radiators  that 
have  more  than  one  feed  port,  such  as  dual-polarized  and  multiple-frequency  antennas.  All  three 
codes  would  benefit  from  faster  matrix  solution,  which  could  result  from  improved  iterative 
methods,  perhaps  using  preconditioning. 

The  third  and  final  category  of  recommendations  deals  with  related  or  similar  electromag- 
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netic  problems  that  may  benefit  from  application  of  hybrid  finite  element  methods.  One  that  is 
a  particularly  straightforward  extension  of  the  present  work  is  gratings  and  frequency  selective 
surfaces.  That  problem  may  be  addressed  simply  by  replacing  the  waveguide  feed  with  a  second 
periodic  radiation  condition  and  an  incident  field  in  the  form  of  a  plane  wave.  A  more  difficult 
extension,  but  one  with  ruany  practical  applications,  is  a  finite-by-infinite  array.  Such  a  model 
could  be  used  to  predict  the  radiating  properties  of  line-source  arrays  or  to  assess  edge  effects  in 
planar  arrays.  Last,  the  calculation  of  scattering  from  objects  that  include  cavities,  or  of  coupling 
into  circuitry  inside  those  cavities  are  problems  that  may  benefit  from  the  use  of  the  finite  element 
method  to  model  the  cavity  interior,  with  an  integral  equation  boundary  condition  at  the  aperture. 
These  problems  are  more  easily  addressed  now  because  the  present  work  has  provided,  among 
other  things,  a  body  of  well-tested  computer  routines  for  three-dimensional  finite  element  and 
waveguide  mode  computations. 
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Appendix  A 

The  Electric  Field  Functional 


A.l.  Variational  Principle  vs.  Weak  Form 

Publications  dealing  with  electromagnetic  finite  element  applications  usually  begin  with 
one  of  two  forms  of  functionals1  with  little  or  no  discussion  as  to  why  one  is  used  and  not  the 
other.  This  appendix  discusses  their  origins  and  the  circumstances  in  which  each  is  appropriate. 
The  conclusion  is  that  the  results  of  the  two  methods  are  indistinguishable  for  typical  electromag¬ 
netic  radiation  and  scattering  problems. 

The  first  alternative  is 

n  l  r 

(see,  for  example,  Jin  &  Volakis  [53])  where  fi  is  the  volume  region  over  which  the  unknown 
field  is  to  be  found  and  T  is  its  enclosing  boundary.  The  subscript  v  denotes  that  this  functional 
is  the  variational  principle ,3  The  second  alternative  is 

FW(E)  =  JjJ  JLvx  W -VxE-k20erW -E 
0  IA 

(see,  for  example,  D’Angelo  &  Mayergoyz  [24]).  Here,  W  is  a  trial  function  whose  form  is  yet 
to  be  determined.  The  subscript  w  on  F  denotes  the  fact  that  this  is  the  weak  form. 

2  A  functional  is  a  mapping  from  a  space  of  functions  to  the  complex  numbers.  It  is  usually 
an  integral  containing  an  unknown  function.  The  result  of  the  integration  is  a  single  number,  in 
contrast  to  an  integral  equation,  which  maps  the  function  into  another  function. 

3  A  variational  statement  may  be  either  a  variational  principle  or  a  weak  form. 


97 


The  variational  principle  Fv  has  been  deduced  from  a  general  from  of  energy  functional 
[54]  and  is  used  in  conjunction  with  the  Rayleigh-Ritz  principle  in  order  to  form  a  system  of 
equations.  The  weak  form  is  derived  more  directly  by  simply  taking  the  inner  product  of  the 
operator  equation  with  the  trial  function.  It  will  be  used  in  conjunction  with  the  method  of 
weighted  residuals  to  form  the  system  of  equations.  In  the  case  of  Galerkin’s  method  (a  special¬ 
ization  of  weighted  residuals)  the  two  forms  may  give  exactly  the  same  system  of  equations.  But 
in  order  to  use  Galerkin’s  method,  the  expansion  function  that  is  admissible  in  the  original 
problem  must  also  be  admissible  in  the  adjoint  problem. 

The  following  section  will  discuss  the  meaning  of  the  adjoint  problem  and  what  the 
conditions  are  for  self-adjointness  as  applied  to  the  vector  wave  equation.  The  third  will  then 
show  that  typical  waveguide  continuity  conditions  represent  non-self-adjoint  boundary  conditions. 
The  last  section  will  show  that  under  the  assumption  that  Galerkin’s  method  is  applicable,  the  two 
formulations  generate  identical  systems  of  equations. 

A. 2.  The  Adjoint  Problem 

Before  deciding  whether  to  use  (A.  1)  or  (A. 2)  one  must  know  whether  or  not  the  operator 
equation  is  self-adjoint.  If  not,  then  the  properties  of  the  adjoint  problem  must  be  determined  in 
order  to  ensure  that  the  trial  functions  are  capable  of  representing  its  solutions. 

The  operator  equation  is  the  time-harmonic,  source-free  vector  wave  equation  for  the 
electric  field  in  a  linear,  isotropic,  inhomogeneous  region: 

L{E)  =  Vxi  VxE-kltl  =  0  (A. 3) 

Hr 

Its  inner  product  with  an  arbitrary  complex  function  W  is 
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<L(E),W > 


VxE-klerE 


•  W*  dv  =  0 


(A.  4) 


Using  a  Green’s  identity  (integration  by  parts)  twice  shifts  both  derivatives  from  E  to  W: 


<L(E),wy  = 


If  [t, 

0  L 


VxE-  VxW* 


~ k\trE . 


-  J|  J-  W*  x  V  x£  • 
r 


(A.  5) 


VxJ_VxW-iJc*  W 


c/v 


Ilf- 

0 

||  _L[w*  xVx£-£xVxiv*]-  fids 


(A. 6) 


From  the  definition 


<£(£),  W>  =  <£,£“(  W)>  (A/7) 

it  is  evident  that  the  term  in  brackets  in  the  first  integral  of  (A. 6)  is  L\  the  adjoint  operator.  It 
is  simply  the  wave  equation  with  the  constitutive  parameters  replaced  by  their  complex  conju¬ 
gates.  It  is  now  clear  that  W  must  be  in  the  domain  of  the  adjoint  operator. 

The  definition  of  self-adjointness  is  L  =  La,  which  obviously  cannot  be  true  if  the  problem 
includes  lossy  materials.  But  it  also  depends  on  whether  the  boundary  conditions  are  such  as  to 
make  the  surface  integral  in  (A.6)  vanish,  i.e. 


|J  2.  Ea  *  x  V  x  £  •  fids  =  ||  —  £  x  V  x  £"  *  •  fids 
r  r 


(A.  8) 


* 

|J  2j-Ea*  xH>  Ms  =  ||  2±-ExHa 


fids 


(A. 9) 
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Suppose  that  all  lossy  magnetic  materials  are  confined  inside  Q  so  that  along  the  boundary  T,  the 
permeability  is  entirely  real.  Then  (A.9)  is  satisfied  when  Ea=E*  and  Ha=H*.  Recall  that  the 
time-harmonic  and  time-dependent  fields  are  related  by  [55:15] 

E(x,y,z,t)  =  y/2R  e(Ee^ot)  (A10) 

(boldface  represents  the  time-dependent  quantity;  the  expression  for  magnetic  field  is  identical). 
Evaluating  (A.  10)  at  any  point  r  along  T,  with  E0  denoting  E(r): 


E(r,t )  =  \/^[Re(£0)coswr  -  Im(£0)sinwf]  (A. 11) 

and  assuming  that  Ea=E* 

Ea(r,t )  =  e,u,J  =  ^2TjRe(£0)cos«r  +  Im(£0)sinwr] 

=  V^Re[£0^‘;"']  (A.  12) 

=  E(?,-t) 


In  other  words,  the  adjoint  fields  are  time-reversed  versions  of  the  original  fields.  They  carry 
power  across  the  boundary  in  the  opposite  direction  and  they  encounter  materials  that  have  gain 
instead  of  loss.  This  is  the  physical  interpretation  of  the  adjoint  problem,  and  is  consistent  with 
the  property  that  if  the  operator  L  is  causal,  then  the  operator  La  is  anti-causal  [56:356]. 

Notice  that  if  the  boundary  T  is  comprised  entirely  of  perfect  electric  and  perfect  magnet¬ 
ic  conductors,  then  (A.9)  is  always  satisfied  because  there  cannot  be  any  transfer  of  power  across 
such  boundaries.  This  leads  to  the  suspicion  that  the  boundary  conditions  that  will  cause  non- 
self-adjointness  are  those  open  boundaries  where  conditions  of  field  continuity  are  to  be  enforced. 
Section  A. 4  will  demonstrate  that  this  is  indeed  the  case  for  the  waveguide/cavity  apertures  that 
are  considered  in  the  main  body  of  this  report. 
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A. 3.  Continuity  Conditions  for  Waveguide  Apertures 

For  the  waveguide  continuity  problem,  the  boundary  T  reduces  to  the  waveguide  aperture 
rA.  Taking  Ea=E*  and  Ha=H*,  the  following  form  is  equivalent  to  (A.9): 


||  E*  •( AxH)ds 


||  E  •  (fixH*  )ds 

r, 


(A.  13) 


(again  assuming  that  p.T  is  real  along  T).  Consider  a  waveguide  whose  axis  is  the  z  axis  and 
which  joins  the  volume  fi  through  an  aperture  in  its  end  wall  at  z  =  0.  It  is  assumed  to  be  match- 
terminated  at  z<  <0.  The  aperture  field  due  to  the  dominant  mode  and  its  conjugate  are 


E  =  g0(l^C0) 

E*  =  jfoO  +C0*  ) 


(A.  14) 


H  =  zxg0r0(l-C0) 
w*  *  txg0Y0(l-C*) 


(A.  15) 


where  g0  is  a  transverse  mode  function,  Y0  is  the  modal  admittance  and  C0  is  an  unknown 
coefficient.  The  outward  normal  to  fi  is  ft  =  -2,  so  the  left  and  right  sides  of  (A.  13)  give 


||  £  *  ’( nxH)ds  =  T0(l+c0*)(i-c0)|j  \g0\2ds 

ra  r. 


|J  E-(fixH*)ds  =  T0(l  +C0)(1  -C0*  )  ||  |g0|2^ 


(A.  17) 


These  two  are  only  equal  if  C0  is  real,  which  is  not  generally  the  case.  Therefore,  the  continuity 
condition  across  a  waveguide  aperture  will  render  the  boundary  value  problem  non  self-adjoint; 
the  functional  (A.l)  is  not  appropriate  even  when  there  are  no  lossy  materials;  and  the  weak  form 
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(A. 2)  should  be  chosen. 


A. 4.  Galerfan’s  Method  vs.  Rayleigh-Ritz 

In  the  main  body  of  this  report,  the  finite  element  method  is  used  to  produce  a  system 
of  equations  from  the  weak  form  functional.  The  field  is  represented  by  a  summation  of  un¬ 
known  complex  coefficients  with  known,  linear  vector  functions: 

£  «  £  =  52  es\J/s(x,y,z )  (A-18) 

S  =  1 


For  the  expansion  functions  to  be  admissible,  they  must  be  in  the  domain  of  the  functional. 
Linear  functions  meet  that  requirement  since  their  first  derivatives  are  continuous  (integrable). 

The  residual  error  is  defined  as  R=L(E)-L(E).  N  weighting  functions  will  be  chosen, 
each  of  which  is  orthogonal  to  the  residual  so  that  (R,wr)=0,  giving  the  system  of  N  equations 


+  Jko1lo  ff  (wrxH)  'fids,  r-  \  ,2,...,N 
r 


(A.  19) 


The  choice  wr=$r  satisfies  the  orthogonality  requirement,  but  it  is  also  required  that  \pT  be  an 
admissible  expansion  function  for  the  adjoint  electric  field.  Starting  with  the  adjoint  operator 
equation,  forming  (E,La(W))  and  using  the  Green’s  identity  once  gives 


<E,La(W )>  = 


V  X  E  -  &q€ 


fjj  [ivx*’ 

-  JJ  J_£x  Vx  W*  •  fids 


rW*  •  E 


dv 


(A.20) 


which  makes  it  evident  that  expansion  functions  for  W  must  have  continuous  first  derivatives. 
Therefore,  the  same  expansion  functions  are  admissible  in  both  the  original  and  the  adjoint 
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problems.  It  is  this  fact  that  allows  Galerkin’s  method  to  be  employed. 


Consider  the  variational  functional  (A.l)  with  the  expansion  (A.  18)  for  the  electric  field: 


1(1 


r= 1  *=1 


fl  L 


±Vx$r- Vxts-k20ertr-$s 
Pr 


dv 


N 

+  jkoVo'52  er  ||  ($rxH)'Ads,  r=l,2,...,Af 

r=l  fJ 


(A. 21) 


The  Rayleigh-Ritz  method  equates  the  stationary  point  of  Fv  to  the  minimization  of  the  above 
with  respect  to  each  of  the  coefficients  er,  i.e. 

dF (£) 

8FV(E )  =  0  «•  — v—~  =  0  ,  Vr  (A.22) 

3er 


Carrying  out  the  partial  derivatives  in  (A. 21)  gives 


N 

»  -  E*. 

s= 1 


tt  L 


J_Vx$r*  VY.$s-k2Qeryr  - 


Mr 


dv 


JkoVo  \\  (JrxH)-flds,  r  =  \  ,2, ...  ,N 


(A. 23) 


which  is  identical  to  (A.  19)  with  wr  replaced  by  yj/T.  This  shows  that  the  systems  of  equations 
resulting  from  the  variational  principle  and  from  the  weak  form  are  identical.  Thus,  under  at 
least  some  circumstances  the  distinction  between  the  two  forms  is  inconsequential,  and  the 
variational  principle  may  also  be  used. 
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Appendix  B 

Waveguide  Mode  Function  Inner  Products 

B.l,  Approach 

The  waveguide  interaction  terms  require  the  computation  of  two  surface  integrals  $si  and 
¥si  from  (36)  and  (37)  where  s  and  i  are  the  edge  and  mode  indices.  Each  may  be  found  by 
summing  the  contributions  from  the  individual  triangles  that  share  edge  s,  e.g.  for  triangle  k: 

At 


Ls  is  the  edge  length,  ft  is  the  outward  surface  normal  (into  the  waveguide),  Ak  is  the  triangle 
area,  fj  and  f2  are  the  linear  scalar  finite  elements  associated  with  the  nodes  bounding  the  edge 
and  T  is  the  3x3  simplex  transformation  matrix  for  the  triangle.  When  the  integral  is  trans¬ 
formed,  dx  dy  -*  2Akdtidt2  and 

*;*>  -  L,(*  ■n[-T2)G„.Tl}Gsi*TuO„-T,2Gy2}  <®'2> 


1  1  -,1 

Gvj  =  J dtx  J  giv  tj  dt2  ,  v  =x,y ;  j  =  1 , 2  (B.3) 

o  o 


Similarly, 

*«?  -  L,[T210I,i-TllGa*T1,G„-T,3G,2]  ®-4) 

The  generic  procedure  used  to  compute  the  terms  of  the  matrix  SJE  is  outlined  in  Figure  B1 .  The 
bulk  of  the  computation  is  in  step  2.b.ii,  where  the  integrals  G„j  are  calculated.  In  the  case  of 
rectangular  waveguide,  they  may  be  evaluated  in  closed  form.  For  the  other  two  waveguide 
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FOR  EACH  APERTURE  (A,B): 

FOR  EACH  MODE  (i): 

1.  COMPUTE  PROPAGATION  CONSTANT,  MODAL  ADMITTANCE 

2.  FOR  EACH  TRIANGLE  (k)  IN  APERTURE: 

a.  COMPUTE  SIMPLEX  TRANSFORMATION 

b.  FOR  EACH  EDGE  (s)  BORDERING  THE  TRIANGLE: 

i.  DETERMINE  EDGE  VECTOR  ORIENTATION 

ii.  COMPUTE  G^  FOR  p=x,y  and  j  =  l,2 

iii.  ADD  CONTRIBUTIONS  and  *si(k>  to  4>si  and  *si 

3.  FOR  EACH  EDGE  (s)  IN  APERTURE: 

a.  FOR  EACH  EDGE  (t)  IN  APERTURE: 
i.  ADD  [  j  kfl  7j0  Y;  <f>si  ^si  ]  to  Sst,E 


Figure  Bl.  Procedure  for  Matrix  Fill  Calculations  Involving  Waveguide  Modes 


types,  circular  and  circular  coaxial,  they  are  computed  numerically  using  Gaussian  quadrature. 
The  remainder  of  this  appendix  discusses  the  details  of  those  integrations.  Expressions  for  the 
mode  functions,  cutoff  wavenumbers  and  modal  admittances  may  be  found  in  [17]  and  [28]. 

B.2.  Rectangular  Waveguide 

The  mode  functions  for  rectangular  waveguide  will  have  indices  m,n  and  p,  where  p- 1 
or  2  for  TE  or  TM,  respectively.  The  vector  components  are 

-  Cmnpx  cos^sin^)  (B.5) 

sin(=If)cos(^)  (B.6) 

where  a  and  b  are  the  waveguide  dimensions  along  the  x  and  y  axes,  respectively.  The  normal- 
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iztion  coefficients  ensure  that  the  modes  are  orthonormal  over  the  waveguide  cross  section. 

The  inverse  of  the  simplex  transform  matrix,  T1,  will  give  x  and  y  in  terms  of  t,,  t2  and 
constant  coefficients: 


jc  =  a0  +  a,  tx  +  a2  f2  (B.7) 

y  =  0O  +0,/,  +jS2/2  (B.8) 

Let  and  represent  the  following  combinations  of  and 

t  _  W7ra(  ^  (B.9) 

„  _  m-irctj  nirfij  (B  W) 

Vi  a  b 

In  terms  of  these,  the  mode  functions  may  be  rewritten: 


=  jsin(^0  +  ^1r)+^2r2)-Isin(T]0  +  7j1r1+Tj2r2)  (B.ll) 

mnpx 

plUEl  =  isin(^0  +  ^1r1  +  ?2r2)  +  isin(r?0  +  r?1r1+?/2r2)  (B12) 

^ mnpy 

Let  H|  denote  the  integral 

i  i-'i 

//^  =  i  |  dtx  |  sin (B13) 
o  o 

Hjj  is  identical  with  iji  replacing  ^  everywhere.  H^’  and  H^’  are  the  same,  but  with  £,  and  £2 
or  77,  and  ij2  reversed.  In  terms  of  these,  the  integrals  in  (B.3)  are 


0,1  -c„,AHrH,) 


(B.  14) 


Gy\  =  Cmnpy{Hi  +  H  q) 

(B.  15) 

<3,2  -  C.,„(hJ-h') 

(B.  16) 

0,2  -  c„w(h'.h',) 

(B.  17) 

Finally  may  be  evaluated  in  closed  form,  but  it  must  be  accomplished  separetely  for  several 

special  cases,  as  given  below: 

(a)  general  case: 

HS  =  -^[cos(*o  +  *,)-  costo^jsinUo^,)] 

2  2i 

(B.  18) 

- - _[cos(£0  +  £,)  -eos(£0  +  £2)  +  (£,  -£2)sin(£0  +  £i)]J 

J 

(b)  £2;*0,  Sj^O,  | 1 1  **1  (for  small  values  of  £2-£j  (B.  1 8)  is  numerically  unstable, 
susceptible  to  large  roundoff  error  and  overflow): 

**1  *  2^-  l^i  2[cos(^o ^ i>  ~  cos$0  +  £,sin($0  +  $,)] 

-  j  cos(|0  +  $2)  +  i({,  -  £2)sin(£0  +  |,)j 

(c) 

Hk  =  -^4{(4^+6{2-fr^)cos(io)-(4?I^6f2)cos(i0^1) 

4*i 

~  (2$f +4£,£2)sin(£0)  -  2(£^  +  £,  £2)sin(£0  +  £, )} 


(B.  19) 


(B.20) 
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(d)  IfiNl: 


Ht  «  —Lj  icos^o-i^sinSo 

2€2  l 


- - —  [ c°s ( %q  +  $,)  -  cos(£0  +  $2)  +  (£]  -  £2)  sin(£0  +  ii) ] 

(e)  \i2\<\,  |^,|  <  1: 


(B21) 


Hk  =  -L|sin€0  +  (i  +  ^)cos«0)J  &21> 

Note  that  the  last  four  are  correct  in  the  limits  as  those  quantities  specified  as  much  less  than 
unity  go  to  zero.  Precedence  is  given  to  the  five  formulas  in  reverse  order,  checking  for  condi¬ 
tion  (e)  first,  and  executing  (a)  only  when  none  of  the  others  (d),  (c)  or  (b)  are  true. 


B.3.  Gaussian  Quadrature  Integration 

For  circular  and  coal'd  mode  functions  the  inner  products  (B.3)  may  not  be  accom¬ 
plished  in  closed  form,  so  they  have  been  evaluated  numerically  using  Gaussian  quadrature. 
Quadrature  formulas  only  apply  strictly  to  one  dimension,  but  are  easily  extended  to  two  dimen¬ 
sional  integration  over  rectangular  areas.  To  use  these,  the  triangle’s  geometry  is  transformed 
to  a  unit  square,  using  an  approach  suggested  by  Stroud  &  Secrest  [57],  The  transformation  to 
simplex  coordinates  mapped  an  arbitrary  triangle  into  one  with  vertex  coordinates  are  (0.°V  (0,1) 
and  (1,0)  in  (tj,t2)  coordinates.  The  second  transformation  is  given  by 

'2  t  -zf  1 

1  (B.23) 

1  ,  t,  =  1 

The  Jacobian  of  this  transformation  is  (1-tj)1,  so  that  a  typical  integral  term  transforms  into: 
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(B.24) 


l  *~'i  l  l 

G,  =  jt1dti  J  g(/,,/2)</t2  =  Ju,dM  J  (1  -u2)2g(ul,u2)du2 
oo  oo 

Let  uk  and  um  denote  the  one-dimensional  quadrature  sample  points  along  ut  and  u2,  respectively, 
with  wk  and  wm  the  corresponding  weights.  Then  the  integral  is  approximated  by  the  sum 

Q  Q 

G\  *  £  Wk“kY,  Wm^-Um)2  S{Uk,Um)  (B-25> 

i=l  m=\ 

where  Q  is  the  order  of  the  quadrature  formula.  See,  for  example,  [57]  or  [58:887]  for  tables 
of  weights  and  quadrature  points. 
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Appendix  C 

The  Periodic  Integral  Equation 

Conventional  derivations  for  the  scanning  properties  of  phased  array  antennas  are  often 
given  in  terms  of  "Floquet  modes,"  which  are  essentially  plane  waves  propagating  in  several 
discrete  directions  away  from  or  towards  the  array.  That  derivation  is  an  analogue  to  ordinary 
waveguide  modes,  and  was  developed  as  convenient  means  for  deriving  mode-matching  solutions 
to  radiation  from  waveguide  arrays.  This  appendix  presents  an  alternative  derivation  using 
Fourier  transforms  that  may,  in  some  cases,  lead  to  greater  insight  than  the  mode-matching 
solution. 

This  derivation  proceeds  from  the  integral  equation  for  an  arbitrary  sheet  current.  It  is 
then  specialized  to  the  case  of  an  infinite  periodic  sheet  current  through  the  periodicity  condition. 
Its  spatial  frequency  domain  representation  is  then  obtained  through  straightforward  application 
of  Fourier  transform  theorems  and  a  Fourier-Bessel  transform  for  the  free  space  Green’s  function. 
The  inverse  transform  then  yields  the  desired  result,  which  is  an  expression  for  the  integral 
equation  not  as  a  continuous  integral,  but  as  a  summation  over  sample  points  in  spatial  frequency. 
An  extension  of  this  derivation  for  skewed  array  lattices  is  also  provided.  The  results  of  Galer- 
kin’s  method  are  specialized  to  the  linear  vector  finite  element  functions  and  analytic  expressions 
for  the  resulting  inner  products  are  given.  Finally,  the  derivation  in  terms  of  Floquet  modes  is 
shown  to  provide  an  identical  system  of  equations. 

C.  7.  The  MFIE  for  Planar  Current  Sources 

The  objective  of  this  section  is  to  obtain  an  integral  equation  for  the  fields  in  the  half 
space  above  the  radiation  boundary  due  to  the  fields  below  it.  The  use  of  the  equivalence  princi¬ 
ple  will  simplify  the  derivation.  The  boundary,  which  will  be  taken  to  be  located  at  z=0. 
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supports  equivalent  currents  M  and  J.  The  source  of  these  are  the  tangential  fields  just  below 
the  boundary  : 


M  =  ixE{z  =  0.) 

(C.l) 

J  =  -z  xH(z= 0_) 

(C.2) 

The  equivalent  problem  in  the  z>0  half  space  also  sees  a  conducting  boundary,  but  it  supports 
-M  and  -J.  This  equivalent  current  is  the  source  of  an  electric  vector  potential,  F: 


oo 

F(7)  «  -±-  j  j  -M(r)G(r-r' )dx' dy'  (C.3) 

-OO 


00 

A(~r)  =  £  ^-J(j)G(?-7')dx'dy'  (C-4) 

-Oo 


G(r~r' ) 


(C.5) 


where?'  and  r  denote,  respectively,  source  and  observer  coordinates,  and  G  is  the  time-harmonic 
free  space  Green’s  function.  Here,  r'  is  confined  to  z=0,  but  r  may  be  anywhere  above  z=0. 
The  magnetic  field  at  the  observation  point  is  [59:36]: 

H( 7  )  =  -JuF  +  -—  I  (C.6) 

jape 

An  integral  equation  is  obtained  by  applying  a  boundary  condition  to  the  above  radiation  integral. 
Specifically,  the  total  tangential  magnetic  field  at  z=0  is  H  =  lx  J: 


flxj  = 


1 


J 


k2F  * 


d2Fx  d2  F 


dxl  dxdy 


k~Fy  + 


d  2F„ 


d2Fr 


dy~  dxdy 


(C.7) 
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Note  that  an  electric  current  source  in  the  z=0  plane  produces  a  magnetic  field  that  is  entirely 
normal  to  that  plane.  Hence  A(z=0)  =  2AZ  ,  so  J  does  not  contribute  to  Ht(z=0)  and  (7)  is 
a  complete  expression  for  the  integral  equation.  The  next  section  will  specialize  this  to  the  case 
in  which  M  is  an  infinite  periodic  source. 

C.2.  The  Periodic  Magnetic  Field  Integral  Equation 

According  to  Floquet’s  theorem,  the  fields  and  currents  anywhere  on  and  above  an 
infinite  periodic  array  must  obey  the  relationship 

$(x+mdx,y+ndy)  =  $(jc,y)e  ^yndy  (C-8) 

where  $  may  be  any  of  E,H,M,  or  F,  dx  and  dy  are  the  lattic  spacings  in  x  and  y  and  \j/x  and  ^y 
are  the  phase  shifts  necessary  to  produce  a  plane  wave  propagating  in  the  %</> 0  direction: 

\{/x  =  k  sin cos^o  (C.9) 

\py  =  k  sin0o  sin0o  (C.10) 

Let  Muc(x,y)  denote  a  unit  cell  magnetic  current  that  is  equal  to  the  source  distribution  within  the 
rectangular  area  -dx^x<  dx  and  -dy  ^y  <  dy  and  zero  elsewhere. 

Consider  an  infinite  two-dimensional  sequence  of  Dirac  delta  functions  located  at  lattice 
points  x=mdx,  y=ndy  for  -oo  <m,n<  ».  Figure  Cl  illustrates  that  the  effect  of  a  two-dimen¬ 
sional  convolution  of  this  Dirac  sequency  with  the  unit  ceil  current  distribution  is  to  replicate  the 
current  distribution  around  each  of  the  lattice  points.  If  each  impulse  is  also  weighted  by  the 
complex  exponential  representing  the  beam  steering  phase,  then  M  is: 

M(x,y)  =  MUC*Y,  {x-tndx,y-ndy)e  j**x e  (C.  1 1) 

m  n 

where  *  denotes  the  convolution  operation.  This  is  an  alternative  form  of  Floquet’s  theorem. 
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UNIT  CELL 
FIELD 


TOTAL  FIELD 


Figure  Cl.  Two  Dimensional  Convolution  of  Unit  Cell  Field  with  Dirac  Impulse  Sequence 

(All  summations  may  be  assumed  to  have  limits  of  -oo  to  +  <»  unless  otherwise  noted.) 

Let  M  denote  the  two-dimensional  Fourier  transform  of  the  magnetic  field  with  respect 
to  the  spatial  frequency  coordinates  kx  and  ky  (underbar  will  be  used  to  denote  transformed 
quantities): 


M(kx,ky)  =  J  J  M{x,y)eiKx  ejk>y  dxdy 


(C.12) 


M(x,y)  =  -L  J  J  M(kx,ky)e ~ik*x e~jk»y dkxdky 


(C.13) 


The  following  four  properties  of  Fourier  transforms  are  required  (see,  for  example  [30: 199-200]): 

&{f*g}=fg  (c-14> 
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&{F(x)e  -jux)  =  F(kx*  a) 


(CAS) 


W-] 

1 dx  J 

•  =  ~J\E 

(C.16) 

H  E 

5  ( x-nd )  } 

-T  £  6(*^) 

(C.17) 

(The  last  one  is  a  form  of  the  Poisson  sum  formula.)  Substituting  (C.9)  into  (C.10)  and  using 
properties  (C.14),  (C.15)  and  (C.17)  results  in 

mx,ky)  =Mc^EE^-4,vy  (ci8> 

x  >  m  n 


2-irm 


-'Px 


(C.19) 


k  =  111-* 
Kyn  Yy 

y 


(C.20) 


The  points  kxm  and  kyn  are  sample  points  in  the  spatial  frequency,  or  spectral  domain.  They  may 
also  be  recognized  as  the  so-called  Floquet  harmonics. 

When  the  source  and  observer  arc  both  in  the  z  =  0  plane  the  Green’s  function  is  only  a 
function  of  x  and  y  and  the  integral  (C.3)  is  a  convolution  integral  written  as 


F(x,y)  =  -±-  M(x,y )  *  G(x,y) 


(C.21) 


whose  Fourier  transform  is 


F(kx,ky)  =  --LMG(kx,ky)  (C.22) 

The  transform  of  the  Green’s  function  is  found  using  a  Fourier-Bessel  transform  [60:12],  with 
the  remarkably  simple  result: 
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or  in  terms  of  the  unit  cell  electric  field: 

i-  EE  T  ■luMk>-k^’k,~k,n) 

m  n 


(C.26) 


(C.27) 


~  =  1  ^  fy  k*ky  (C.28) 

2k K'  [  kxky  (* 2-k2x) 

Finally,  testing  is  to  be  carried  out  in  the  spatial  domain,  requiring  the  inverse  transform  of 
(C.25).  With  the  delta  function  in  the  integrand,  that  integration  reduces  to  a  sampling  of  the 
integrand  at  each  and  kyn: 

j  -  E£  Tmn  ■  EJkxm,kyn)e -Jl-Xe (C. 29) 

m  n 

This  form  of  the  MFIE  is  only  valid  for  rectangular  array  lattices,  but  many  actual  phased  arrays 
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m”  2  dxdykKmnr, 


tf'fyn)  kxmkyn 
kxmkyn  (k2-klj 


(C.30) 


k 2  -  k2xm  -  k2 

xtn  yn 


1/2 


use  triangular  lattices.  The  following  section  will  show  how  (C.29)  is  modified  to  accomodate 
skewed  lattices. 


C.3.  Skewed  Array  Lattices 

Phased  array  antenna  elements  are  usually  arranged  in  a  triangular  lattice,  formed  by 
shifting  successive  rows  to  the  right  or  left  by  one  half  the  column  spacing.  This  allows  a  larger 
inter-element  spacing  (hence  fewer  elements  for  a  given  aperture  area)  to  cover  a  given  scan 
region  without  grating  lobes. 

Figure  12  shows  an  even  more  general  case  in  which  the  shift  between  successive  rows 
is  not  necessarily  one  half  the  column  spacing  dx.  The  Floquet  condition  for  this  situation  is 

£(jc +mdx  +/i dycot7  ,y  +  ndy)  =  £(jt,y)e  j'f'*imd**ndycolt)e  j%ndy  (C.31) 
or  in  convolution  form 

E(x,y)  =  Euc(x,y)  *  £  £  8(x-mdx-ndycoty ,y-ndy)e  j^xe  j*>y  (c.32) 

m  n 

The  Fourier  transform  of  (C.30)  is  required,  but  it  may  be  obtained  without  directly  performing 
the  integration.  As  in  the  case  of  the  rectangular  lattice,  a  result  similar  to  (C.  16)  is  expected, 
i.e.  the  unit  cell  transform  times  a  series  of  spectral  domain  delta  functions.  The  Fourier  trans¬ 
form  pair  (C.9)  and  (C.  16)  have  a  unique  interpretation  in  terms  of  direct  and  reciprocal  lattices 
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[61:94-98].  The  coordinates  (kx,ky)=(2irm/dx,2irn/dy)  are  points  in  the  reciprocal  lattice  corre¬ 
sponding  to  (x,y)=(mdx,ndy)  and  (4ir2/dxdy)  is  the  area  of  one  of  its  unit  cells. 

In  the  skewed  lattice,  the  coordinates  of  any  element  are  integer  multiples  of  the  basis 
vectors  a  and  b: 


a  =  idx 

b  =  Jtd  coty  +  yd 


(C.33) 


The  basis  vectors  in  the  reciprocal  lattice,  a  and  0  are  found  by  solving 


with  the  result 


a  •  a  -  b  •  i 3  =  1 

a  •  /8  =  b  •  a  =0 


(C.34) 


a  =  J? 


1 

coty 

dx 

J 

d* 

(C-35) 


0  =9 


(C.36) 


The  spatial  frequency  coordinates  corresponding  to  points  in  the  reciprocal  lattice  are 


k'xmn  =  '(ma+nP)  = 


(C.37) 


k'ymn  =  2irj>-(ma  +  /I/J)  =  111  -  2™Sot7  {C.38) 

ay  ax 

The  primes  signify  that  these  are  the  unscanned  lattice  points.  The  unit  cell  area  in  a  skewed 
lattice  is  the  same  as  in  a  non-skewed  lattice  with  the  same  dx  and  dy,  so  the  spectral  domain  unit 
cell  area  is  also  the  same.  The  end  result  is  that  the  Fourier  transform  of  (C.32)  is 
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(C.  39) 


E(kx,kv )  =  E 

—  v  x  >  y/  —  ac 


47T2 


x  y  m 


£  ) 
I'ymri' 


~^Osin^Ocos^O  (C-40) 

kymn  =  ^  "  ---^COt-  ~  *osin0osin<£o  (C.41) 

Uy 

The  effect  of  beam  steering  is  to  shift  all  of  the  lattice  points.  The  last  three  equations  above 
explain  the  origin  of  the  Floquet  harmonics  for  a  skewed  lattice,  and  agrees  with  the  result  given 
by  Mittra  et.  al.  [62:1596].  The  relationship  between  2xH  and  E  is  still  as  given  by  (C.29), 
except  that  the  sample  points  are  now  k^^ky^  instead  of  kxm,kyn. 


C.4.  Expansion  Function  Fourier  Transform 

The  electric  field  within  the  unit  cell  is  expanded  in  known  vector  functions  \j/s  with 
unknown  complex  scalar  coefficients  es.  By  linearity  of  Fourier  transforms,  the  unit  cell  field 
is 


Euc  =  E  «.*,<**.*,)  (C42) 

s  =  1 

where  £8  is  the  two  dimensional  Fourier  transform  of  \JS.  These  may  be  evaluated  analytically 
with  the  help  of  homogeneous  coordinates  within  the  triangles  subdividing  the  radiation  boundary 
(the  faces  of  those  tetrahedra  bordering  on  the  boundary). 

Suppose  that  within  triangle  k  edge  s  goes  between  local  nodes  i  and  j.  Then  in  terms 
of  the  2D  homogeneous  coordinates. 


=  Lijit^tj-tjVt ,) 


(C.43) 
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Let  Tj  denote  the  Fourier  transform  of  the  scalar  linear  finite  element  defined  at  node  i,  denoted 
fj.  In  the  homogeneous  coordinates,  f;=tj  and 


t4  =  J  J  ejk*xejk>y  dxdy  =  2  A  J  dt2  J  tt  ejkx* ejk*y dtx  (C.44) 

-oo  0  0 

The  factor  2A  is  the  inverse  of  the  Jacobian  of  the  transform  from  (x,y)  to  (tj,^.  The  inverse 
coordinate  transform  is  expressed  in  terms  of  six  constants  given  in  Appendix  B,  (B.7),(B.8). 
Let  Bj  denote  the  following  combinations: 

=  aikx  +  fl,ky  ,  /  =  1,2,3  (C.45) 

Then  substituting  into  (C.44): 

1  *  _,2 

r,.  =  2AeiB°  J  ejB2hdt2  j  dtx  (C-46) 

0  0 


Five  separate  cases  must  be  considered,  depending  on  the  range  of  the  constants  B,. 


(a)  general  case 


t j  =  2  AejB° 


~j 

b]b2 


JejB2  +  eiB\jB2-2jBl  + 
B2{BX-B2f  b]{Bx-B2)2 


(C.47) 


(b)  |B,-B2|<<1 


Tl  = 


2AejB° 

B\ 


1  jB 2  ^  jB\  ^  Bj(B\  B2) 

1 

2  2  3  12 

\-eiB'{\ 


~jB\ ) 


(C.48) 
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(c)  |B2(<<1 


-ejB'\j{2Bx  +  3 B2)+B,(B1+B2)] 
fli 

7(2fi,  +  3B2-ifi?B2) -£,(£,  +  2B2)]  J 


(d)  |  B,  |  <  <  1 


T/  =  ™?Ll{EjB'(jB2+yBx)  + 

b24 

(2BxB2  +  B2  -  ^BxB2)  +j(-2Bx  J 


(e)  | B j  |  <  <  1  and  jB^I  <  <  1 


(C.49) 


(C.50) 


r;  =  ±^L_  (20  +  lOyB,  +  5 ]B2  -  2 BXB2 ) 


(C.51) 


The  small  argument  forms  ensure  numerical  stability.  The  expressions  for  r2  are  obtained  by 
reversing  Bi  and  B^ 


C.5.  Integral  Equation  from  Floquet  Modes 

The  following  is  a  more  conventional  derivation  of  the  periodic  integral  equation  using 
"Floquet  modes."  It  follows  the  same  general  procedure  used  in  deriving  the  integral  equations 


for  waveguides. 


The  orthonormal  vector  Floquet  modes  [29:41-42]  are  analogous  to  waveguide  mode 


functions: 


1  [  *kyn  ~Skxm~\  J (*,,„* +* 


vn?)  /TP  1 


&2t 


(djy 


*kxm  +  Sk 


yn 


/ 


,2  ,2 

xm  yn 


ei(kxn,x*kyny >  (TM) 


(C.53) 


The  corresponding  modal  admittances  are 


pmn 


KmnlkOV0  P=1  (TE> 
kOlKmnVo  P=2  (TM) 


(C  54) 


Using  the  waveguide  integral  equation  from  Chapter  IV  (29)  but  expanding  the  mode  sum  over 
three  indices: 


E  EE 

m=0  n=l  p= 1 


Y  Q 

pmn opmn 


8 pmn  ds  -  J  =  0 


(C.55) 


Et  is  the  transverse  (to  z)  unit  cell  electric  field  on  rR.  Due  to  the  form  of  the  complex  exponen¬ 
tial  factors  in  (C.52)  and  (C.53),  the  integral  (C.55)  results  in  factors  involving  its  Fourier 
transform,  1^.  The  TE  and  TM  mode  terms  from  (C.55)  are 


v  t;  (p  n  _  *mn'k0Vo 

^1  mn  8  1  mn  J  ^ t  *  8  1  mn  ^ — 

rR  dxdy(kxm+kyn) 

~  ^—ykxmkyn  ~SKxkxmkyn  +  9Kykx, 


^2mn  82mn  J  ’  82mn  d$ 


ko^mnVQ 

dxdy(k  2xm+k;„) 


8E  k2+£E  krmkvn+yE  kr,„kv„  +  yE  k2 

— x  x m  — -y  xm  yn  j  x  xtn  yn  j  — -y  yn 


(C.56) 


Summing  the  TE  and  TM  modes: 
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(C.58) 


[*£,(4  +*Z/xmkyn  +SExkxmkyn  +?Ey(k0 -k2xm) 

When  this  is  written  in  dyadic  notation,  it  is  clear  that  the  integral  equations  (B.29)  and  (B.55) 
are  the  same. 
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Appendix  D 


Periodic  Boundary  Conditions 
for  the  Finite  Element  Problem 

This  appendix  develops  the  method  for  applying  periodic  boundary  conditions  to  a  finite 
element  problem  in  one  dimension.  Consider  the  two  functions  f(x)  and  h(x)  shown  in  Figure 


Dl,  related  by  a  linear  operator  equation  Lf=h.  Their  magnitudes  are  periodic,  repeating  on 
each  interval  ,  but  each  interval  has  a  progressive  phase  shift  relative  to  the  next: 

f(x  +  nd)  =  f(x)e^H  t01) 

This  may  be  regarded  as  a  periodic  boundary  condition  for  the  function  on  the  interval  [0,d]. 
For  example  purposes  the  method  of  weighted  residuals  will  be  used  to  produce  a  functional: 

F(f)  =  <Lf, tv>  -  <h,w> 

d  (D.2) 

=  |  [ L(f)w *  -  hw*}dx 
o 

The  functions  f,  h  and  w  will  be  represented  as  sums  of  complex  coefficients  (f; ,  h; ,  and  w;) 
times  scalar  basis  functions  tj(x)  where  i  may  range  from  -«  to  +  « .  N+ 1  of  these  expansion 
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Figure  D2.  Expansion /Weighting  Functions 

functions  are  nonzero  within  the  interval  [0,d]  as  illustrated  in  Figure  D2.  f,  and  h,  are  the 
values  f(0)  and  h(0);  and  fN  and  hN  are  the  values  f(d)  and  h(d).  The  functions  t;(x)  are  not 
necessarily  linear  as  shown,  and  do  not  necessarily  extend  over  subintervals  of  the  same  length. 
However,  those  in  successive  intervals  must  be  replicas  of  each  other  with  ti+N(x)=ti+1(x-d). 
After  expanding  the  functions  in  (D.2),  the  derivative  of  F  with  respect  to  each  wj  gives  an 
infinite  tridiagonal  system  with  matrix  and  right  hand  side  entries 

Oo 

sJfi  =  J  LVtixHtjWdx  (D.3) 

-00 

00 

gj  =  j  h(x)tj{x)dx  (D.4) 

-00 

For  example,  the  equations  pertaining  to  the  nodes  -1,0,1, 2,.. ,,N-:  3  are 
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*0,-1 /-I 

+ 

O 

O 

+ 

*0,1/1 

=  So 

•*1,0/0 

+ 

*1,1 /l 

+ 

s  1,2/2 

=  Si 

*2,1 /l 

+ 

*2,2  fl 

+ 

*2,3/3 

=  S2 

SN- 1  ,iV-2  /v-2 

+ 

*AM,Am/v-I 

+ 

sn-i,n/n 

~  Sv-1 

sN,N-\  fw-l 

+ 

sN,n/n 

+ 

SN,N*  1  /v+  1 

=  Sn 

sn*i,n/n 

+ 

*AM,Am/v+1 

+ 

SN*\,N*2  /v+2 

=  Sn*  1 

The  periodicity  conditions  on  the  discrete  coefficients  f  and  g  are 

fi* n-\  =  fie^  '■>  fi-N*\  =fie  ^  (D.6) 

*i*N- 1  =  Si***  ;  ii-N* i  =  *«■*"■'* 

The  matrix  elements  must  also  satisfy  a  periodicity  condition.  Since  tl+,  and  ti+N  are  identical 
for  all  i,  it  is  clear  that  sj  +  1  i+1  =  si+N  j+N.  We  can  rewrite  the  system  so  that  it  only 
involves  the  unknown  values  of  f  and  known  values  of  g  within  the  interval  [0,d): 


SN-\,N-2fN-2e 

+ 

sN-l,N-\fN-le  J * 

+ 

*AM.A//| 

- 

(a) 

SN,N-\  In-  1 e  ^ 

+ 

*l,l/l 

+ 

*l,2/> 

=  Si 

(b) 

*2,l/l 

+ 

*2,2  f2 

+ 

*2,3/3 

=  S2 

(c) 

sN-\,N-lfN-2 

+ 

sN-\,N-\  /v-1 

+ 

*AMJv/|^* 

=  Sn~\ 

(d) 

*n,am/am 

+ 

*1,|/|^ 

+ 

*1.2  k*** 

=  gleJ* 

(e) 

*2.l/l^^ 

+ 

*2,2  he** 

+ 

*2,3  he^ 

= 

if) 

Multiplying  (b)  and  (c)  by  e^  and  subtracting  from  (e)  and  (f),  repsectively,  eliminates  the  latter 
three  from  the  system.  Similarly,  multiplying  (d)  by  e'^  and  subtracting  from  (a)  eliminates  (a). 
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Continuing  the  process  will  eliminate  all  equations  preceding  (b)  and  succeeding  (d)  leaving  only 
a  reduced  system  of  N  equations: 


'N,N-\fN-\e  ** 

+  *1.1 /l 

+  *1,2  fl 

=  St 

*2,1/1 

+ 

...  lo 
> 

+  *2,3/3 

=  82 

(D.8) 

sN-\,N-lfN-2 

+  *am,am/n-i  +  *AM,Jv/l^ 

=  8n-  1 

Suppose  that  now  the  original  problem  geometry  is  truncated  at  the  boundaries  x=0  and 
x=d,  and  the  two  ends  are  "wrapped"  back  on  each  other,  as  shown  in  Figure  D3.  Now  nodes 


Figure  D3.  "Wrapped"  Domain 

1  and  N  are  the  same  point,  as  are  2  and  N  +  l,  etc.  Tli  n  the  inner  product  in  (D.3)  has  new 
terms  for  (i,j)=(l,N-l),(N-l,l),  and  referring  to  (D.8)  it  is  evident  that 

*i,tv-i  =  sN,N-\e  ^  (D.9) 

f  At- 1,1  =  SN-\ ,Ne^ 

This  system  is  no  longer  tridiagonal  because  the  boundary  terms  have  introduced  new  elements: 
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ru 

sl,2 

0 

0  •  •  • 

f2,l 

s2,2 

s2,3 

0  •  •  • 

0 

s3,2 

•*3,3 

•*3,4  *  *  * 

IS\  = 


0 

0 

0 


0  SN,N-I  e  ^ 

0  0 

0  0 


(D.10) 


0  0  0  0 

Sn-un^  o  0  0 


’  *  *  sN-2,N-2  sN-2,N-l 
*  *  *  sN-\,N-2  sN-\,N-\ 


0 

sN-l,N 


If  the  original  (infinite)  matrix  was  Hermitian  (adjoint),  then  sN  N_,  =  s*N_1N,  and  the  new 
system  is  adjoint  as  well.  This  indicates  that  periodic  boundary  conditions  do  not  necessarily 
cause  an  operator  to  become  non-self-adjoint. 
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